Metaparticles: Computationally engineered nanomaterials with tunable and responsive properties

Nanoparticles (NPs) are versatile platforms in a variety of fields, including biomedical,^1–3 cosmetics,^4,5 material engineering,⁶ and the food industry.⁷ In cosmetics, nanoparticles are used as filters of UV radiation in sunscreens or to enhance the delivery of active ingredients in dermal tissues.^4,5 In the food industry, nanoparticles are used as transport platforms of “nutraceuticals” toward specific tissues and for preservation purposes.^8,9 In the biomedical field, nanosized objects loaded with molecules are commonly used as carriers to deliver their content (e.g., nutrients, drugs, and imaging agents) to a desired target.^10,11 Their diversity in size, functionalization, and available platforms, make nanoparticles excellent candidates for therapeutic activity in precision medicine.^3,12 Medical applications of NPs range from cancer therapy, where nanocarriers are used to actively target cancerous cells^13,14 to immunotherapy¹⁵ and genome editing, where inhalable NPs formulations are used as therapeutic agents in cystic fibrosis.¹⁶ Metal–organic nanocages are gaining momentum in the biomedical context, showing reduced cytotoxicity, good in vivo biodistribution, and anti-cancer activities.^17,18 Nanoparticles are good candidates in the rational design of novel materials with mechanical properties that are inaccessible to ordinary materials, i.e., metamaterials.^19,20 NPs have been used to reprogram the deformation behavior of mechanical metamaterials²¹ and to engineer stretchable materials with a broad range of organic or inorganic functionality²² as well as biomedical applications.²⁰ 

Nanoparticles consist of a core platform functionalized with ligands on the surface. These ligands are designed to enhance and modulate the individual and/or collective properties of NPs in specific environments.²³ For instance, DNA-patched nanoparticles have been synthesized to control the self-assembly of colloidal metamaterials with functionality-explicit architectures.²⁴ In the biomedical field, various nanoparticle cores are designed to fit specific applications.¹ Hence, lipid nanoparticles are frequently used for drug delivery purposes^25,26 because of their high stability, biocompatibility, and specificity in the cellular environment. While also improving drug bioavailability,^27,28 polymeric nanoparticles struggle with poor stability and solubility due to their hydrophobic nature, which often limits their drug delivery abilities.²⁹ Inorganic nanoparticles consisting of metallic cores are abundantly used for imaging purposes.³⁰ Their optical and electrical properties make them suitable for specific delivery conditions, such as thermal-based therapeutics³¹ or, more recently, in vivo diagnostics.^32,33

Computer simulations complement experimental research to provide insight into specific chemical and physical processes, e.g., informed ligand design, property prediction, and mechanisms in biological environments. They aid in the development of biological-compatible nano-structures to guide and inform on processes exceeding experimental resolution.^34–36 At atomistic resolution, molecular dynamics simulations provide insight into the fine details of coated nanoparticles.³⁷ They are powerful tools to suggest and investigate the effects of specific chemical modification on the biological environment, e.g., interactions with membranes,³⁸ behavior in the blood flow,³⁹ or even make predictions that can be subsequently experimentally tested.⁴⁰ Gold nanoparticles are among the most studied candidates due to their size (<10 nm), optimized force-field parameters, and ease in introducing ligand modifications.^41–43 However, atomistic detail comes with a computational cost, particularly when studying accumulations of nanoparticles or biological response, which require access to long time- and lengthscales (e.g., membrane associated processes and diffusion).

Coarse-grained models overcome these difficulties by grouping atoms into single beads, thereby reducing the degrees of freedom of a system and making it computationally more affordable when simulating large systems.^44,45 In particular, nanoparticles designed as single beads decorated with interaction patches to mimic the surface moieties are ideal to explore the influence of size and shape,^46–48 ligand concentration, and distribution^49,50 to investigate cellular uptake and accumulation.⁵¹ These minimal coarse-grained models are suitable candidates to help with gaining insight into generic mechanisms while reaching timescales comparable with the experiments. Most commonly, they use rigid representations of the nanoparticle core,⁵² while in reality the core can be deformable (e.g., can consist of proteins, micelles, and lipids). As such rigid models do not capture the complete nanoparticle morphology and its ability to adapt to environmental conditions and overcome biological barriers.⁵³ Hence, although relevant, the mechanical properties of nanoparticles, which influence their properties at individual and collective level, are often neglected in minimal models. For instance, both experimental⁵⁴ and computational studies^55,56 have shown that soft NPs can enable drug delivery through mechanisms that avoid endocytosis, a process requiring biological membrane curvature and breakage, which slows down absorption and introduces significant complications.^57,58 This process reminds us of internalization viral mechanisms that leverage intrinsic flexibility for more rapid cellular entry.^59,60 In this context, tailoring physical properties of individual nanoparticles to adapt to the environment can aid in the development of multilayer metamaterials.⁶¹ Despite substantial evidence linking mechanical properties to functionality and macroscopic properties of nanoparticles, existing coarse-grained model do not take into consideration the intrinsic flexibility of a nanoparticle.

Here, we develop a tunable coarse-grained model able to capture the intrinsic flexibility of (nano)particles, arising from soft core platforms or anisotropic ligand coverage. We introduce a model, in which a (nano)particle is represented as a set of hard spheres interconnected with flexible bonds that can respond to applied external tension, hence adapting to the environment. Due to their modularity, versatility, and adaptability, we refer to the designed nanomaterial as metaparticles (MPs). Their generic properties and tunability, allows finding a wide range of applications, from nano-scale drug delivery to macro-scale metamaterial development. The topologies of the metaparticles rely on symmetrical structures inspired by the structures of fullerenes (Fig. 1). Our results show that the metaparticles show a stress response similar to that of elastomers. In addition, the response is modulated by the delicate balance between size, topology, and area of applied stress. Finally, the deformation and recovery timescales of large metaparticles depend on the mechanisms of deformability and relaxation.

The model developed here is a generic representation of a particle with responsive and adaptable properties, which we refer to as a MetaParticle (MP). We seek to mimic the mechanical properties of nanoparticles with soft cores, such as nanogels,⁶² polymers,⁶³ and nanoparticles made from dendrimers as building blocks,⁶⁴ micelles,⁶⁵ vesicles,⁶⁶ or biomolecular condensates,⁶⁷ and elastomeric materials with characteristics similar to biological tissues.⁶⁸ To achieve this, we introduce a generalized model of a metaparticle consisting of multiple interconnected beads, symmetrically arranged in topologies inspired by the structures of fullerenes. These are symmetric nanostructures that consist of carbon atoms interconnected by single or double bonds to form closed or semi-closed cage-like molecules. For our model, we consider only the closed topologies. These molecules form an interconnected network of 12 pentagons and a variable number of hexagons or, less common, heptagons,⁶⁹ resulting in hollow structures (e.g., spheres, ellipsoids, and spherocylinders) [Fig. 1(a)]. Topologically, this arrangement translates into a network of vertices joined by three edges. In our metaparticle model, each vertex is represented as a hard spherical bead. To introduce the flexibility element, each bond (or edge) has a spring-like behavior [Fig. 1(b)]. To investigate the mechanical properties of MPs with comparable topologies, yet differing in size and arrangement, we chose ten model systems. They range from small metaparticles with twenty beads in their composition to larger MPs with 240 constituents. The smallest model, MP₂₀, consists only of twelve pentagons, while all the other models have, at least, one hexagon in their topology. To account for variability in the shape and topology, we consider metaparticles with different symmetries (Fig. 1). Furthermore, to characterize systems with the same number of beads but different topologies, we use three model systems with different symmetries for the MPs with 36 vertices.

Throughout this work, we use reduced Lennard-Jones units,^79,80 with ɛ and σ representing the units of energy and length, respectively. Time unit is $τ=σ/m/ε$⁠, where m is the particle mass and the temperature is given in ɛ/k_B (k_B is the Boltzmann constant). Each system was initialized by placing a metaparticle in a cubic box under period boundary conditions (40σ per edge). The template conformations for the MPs were fullerene structures with twenty to 240 vertices⁸¹ [Fig. 1(a)].

The beads in the MP interact via bonded and non-bonded interactions. An individual particle interacts with its neighbors via the elastic FENE potential. We adopted the standard choice of bond strength as introduced in the Kremer–Grest model K = 30(ɛ/σ²).⁸² In addition, we explored the effect of K = 5ɛ/σ² and K = 100ɛ/σ² on the different analyzed properties of the MPs. The maximum extensibility is slightly higher than in the Kremer–Grest model,⁸² i.e., r₀ = 1.8σ to account for a longer adaptability range of the metaparticle. In addition, non-bonded beads can interact via a generic hydrophobic potential described by a Lennard-Jones potential, with an interaction strength ɛ_LJ = 1ɛ, $σ̄=2.5σ$⁠, and if they are within a cutoff distance r_c = 2.5σ from each other. The time step was fixed at Δt = 10⁻⁴τ, where τ is the BD time unit, and the systems were evolved using Brownian dynamics simulations using the LAMMPS simulation package.^73–76 The mobility tensor values correspond to the well-known Stokes drag expressions for spherical particles.⁸³ 

In a typical set of simulations, a metaparticle is first relaxed for 10⁵τ in the absence of restraints and only being subjected to Brownian motion. Next, the MP is mechanically deformed by performing pulling simulations for 10⁵τ. For this, forces were applied on the equilibrated structure. Each force was applied on anti-parallel sets of beads, i.e., antipodal hexagons or pentagons. To avoid shearing effects, each system was realigned with the antipodal surfaces aligned along the direction of applied deformation. In the final step, the applied stress is removed and the metaparticle is simulated in the absence of restraints till each specific MPs is fully relaxed.

First, we evaluated the radii of gyration R_g of the metaparticles as the ensemble average over the positions of all beads in an MP over time [Fig. 2(a)]. We found that R_g increases with the number of beads in the MPs, ranging from 2.01 ± 0.006σ for MP₂₀ to 6.74 ± 0.008σ for MP₂₄₀. The MP₃₆ variants are essentially identical in their radii of gyration, despite the different topologies, i.e., 2.68 ± 0.007σ, suggesting that bead arrangement marginally affects metaparticle size.

Next, we investigated the diffusive behavior of the MPs by analyzing their mean squared displacement as $MSD=⟨R(t0+Δt)−R(t0)2⟩=6DΔt+C$⁠. The brackets ⟨⟩ represent the average over all time origins t₀ within the simulation run, D is the metaparticle translational diffusion coefficient, and $C$ is a constant. In agreement with R_g analysis, we find that the translational diffusivity of the MPs reduces with increasing system size, i.e., the smallest MP in diameter; MP₂₀ diffuses the fastest, while MP₂₄₀ is the slowest [Fig. 2(b)]. The three MP₃₆ metaparticles have comparable diffusion coefficients, showing that topology only marginally affects diffusion. Furthermore, the internal flexibility of the MPs, characterized by different values of the spring constant, only marginally impact their radii of gyration and the diffusion (Fig. S1).

By design, our metaparticles have symmetrical topologies inspired by the structures of closed fullerenes (Fig. 1). Hence, each MP consists of 12 pentagons and a variable number of hexagons, and their arrangement is closely linked to shape and topological properties. Pentagons induce curvature, which is essential to close the metaparticle and form a spherical or quasi-spherical shape. The MPs with less than sixty beads have pentagons that share an edge with other pentagons, with MP₂₀ presenting only pentagons in its structure. Systems with more than sixty beads have pentagons surrounded by hexagons following the isolated pentagon rule (IPR) rule.⁸⁴ 

To gain insight into the topological effects on the behavior of the MPs, we first analyzed their sphericity [the resemblance to a perfect sphere, Fig. 2(c)], asphericity (deviations from a perfect sphere, Fig. S2), and acylindricity (resemblance to a cylinder, Fig. S3). The sphericity was determined from the ratio between the volume of a metaparticle, V_MP, and its surface area A_MP⁸⁵ as $36π⋅V2MPA3MP$ and following the isoperimetric quotient definition.^69,86 V_MP and A_MP were calculated using the alpha-shape method^87,88 as implemented in Ovito.⁸⁹ Values of Ψ close to 1 indicate nearly spherical objects, while decreasing values indicate deviations from a spherical shape with the object becoming nearly flat as Ψ approaches 0.

The results show that the larger metaparticles with icosahedral symmetry, i.e., MP₆₀, MP₁₈₀, and MP₂₄₀, approach the most the spherical shape as this arrangement reduces the strain on the bonds and the angles, while seeking to minimize the surface area-to-volume ratio. MP₁₈₀ has the highest resemblance to a sphere, given the homogeneous distribution of pentagons on its surface [Fig. 2(c)]. Compared to its smaller counterpart, MP₆₀, with the same symmetry, MP₁₈₀ consists of more beads, which allow for more favorable packing, even distribution of strain over the bonds, and smoother curvature, thereby better approximating a sphere. The largest MP of the family, MP₂₄₀, has sphericity slightly lower than MP₁₈₀ despite the higher number of beads and implicitly hexagons in its composition. This is an inherent effect of the sphericity estimation using the isoperimetric quotient definition, which is consistent with previous calculations on fullerenes.⁶⁹ In addition, the metaparticles introduced here are flexible structures subjected to Brownian motion and do not account for electronic motion and chemical stability as in the case of real fullerenes. The MPs with D_6d symmetry, i.e., MP₄₈, MP₇₂, and MP₉₆ have lower sphericity than the large icosahedral MPs. Their sphericity correlates with size [Fig. 2(c)]. Hence, the smaller MP₄₈ has the lowest sphericity due to its oblate arrangement, which has a lower number of hexagons in its composition, leading to higher bond strains compared to MP₇₂ and MP₉₆. The MP₃₆ isoforms have comparable sphericities (0.87, 0.85, and 0.87) independently of their symmetries. Hence symmetry, bead packing, and, implicitly, the distribution of hexagons do not have a major impact on the sphericity of the small metaparticle. MP₂₀, the smallest MP with icosahedral symmetry and the outlier among the MPs has lower sphericity compared to the MP₃₆ isofroms. This is due to its exclusive pentagonal composition, which restricts efficient packing. All in all, topological analysis reveals that metaparticles symmetries affect their sphericity to a larger extent than their size.

To investigate the effects of size and topology on the mechanical properties of the metaparticles, we studied the tensile stress response. Briefly, we deformed the individual structures by applying tensile forces on anti-parallel particle sets, i.e., diametrically opposed pentagons or hexagons. First, each equilibrated MP was recentered in the simulation box and the centers of mass of the set of pulling area were aligned along one of the principal axes to facilitate uniaxial deformation. Then, a series of forces were sequentially applied on diametrically opposed hexagons (i.e., previously aligned sets of particles), except for MP₂₀, which is devoid of hexagons. We applied a total of 38 distinct forces on each metaparticle (see the supplementary material for details), cumulating 38 × 10⁵τ sampling per system. The maximum force applied differs per system and was chosen as the highest stress the MP can withstand prior to bond breakage.

The stress–strain curves have a non-linear response [Fig. 3(a)] and are marginally affected by the stiffness of the springs (Fig. S5). At low strains (<≈5%), the stress response is linear and only marginal mechanical response of the metaparticles is observed, with the thermal noise being dominant [Figs. 3(a) and S4].

At higher strains (>≈5%), the stress response shows an exponential behavior. The bonds start to align in parallel to the direction of the applied stress and the MPs deform. It is notable that the metaparticles heterogeneously resist the applied stress prior to bond breakage, i.e., the maximum applied stress varies across the different MPs. Furthermore, the stress–strain profiles vary with metaparticles size and topology and exhibit a response similar to that of elastomeric materials.

To further quantify the elasticity and response of the developed metaparticles to stress, we evaluated their toughness, T, as the integral of the stress–strain curve [Fig. 3(b)]. The highest toughness is measured by the MPs with D_6d symmetry, with values in the 50–60 ɛ/σ³ interval. Their oblate conformation and the direction of applied stress, i.e., along the short symmetry axis, allow the accommodation of higher mechanical deformation. This may explain the lack of dependence of toughness on metaparticle size. The small MPs (R_g < 3 nm) have comparable toughness independently of topology. For the larger metaparticles, particularly among those with icosahedral symmetry, MP₁₈₀ shows the lowest resistance to stress and MP₂₄₀ the highest, suggesting that size plays a secondary role in defining the amount of energy MPs can store. The response to stress also varies with the location of the applied stress, i.e., the area of applied deformation or the arrangement of the pentagons relative to the pulling hexagon.

To understand this effect, we evaluated the average pulling area (A^p) at the ultimate tensile strength (UTS), i.e., at the maximum applied stress that a system can handle prior to bond breakage [Fig. 3(c)]. Two different regions can be identified depending on A^p. The stress applied to areas of ≈7.8σ as in the case of MP₂₀, MP₇₂, and MP₂₄₀, leads to higher resistance with UTS values over 200 ɛ/σ³. It is notable that the largest and the smallest MPs, i.e., MP₂₄₀ and MP₂₀, respectively, withstand the highest applied stress. More variability is observed among the metaparticles with the stress subjected to areas of almost twice the size. Among them, MP₁₈₀ has the lowest UTS, which may contribute to the reduced toughness of the metaparticle and potentially relate to the propagation of stress along its surface. Hence, our results show a correlation between the average pulling area and the ability of the MPs to respond to the applied stress. This correlation is independent of size and topology.

To understand how the induced deformation propagates along the surfaces of the metaparticles, we evaluated the energy distribution per bond at different stages during the deformation. Briefly, we calculated the average energy per bond as the mean over the unconstrained runs and over the pulling simulations, evaluated at the highest strain values (Fig. 4). The analysis focused on the systems devoid of external forces reveals that generally the overall energy per pentagon is constant, while the energy is distributed differently per bond depending on the neighboring polygon [Fig. 4(a)]. For instance, the energy per pentagon bond in MP₄₈ is higher when the edge is shared with another pentagon, than when it is shared with a hexagon. By comparison, the energy per bond is homogeneously distributed in the systems where a pentagon is surrounded only by hexagons [e.g., Fig. 4(a) for MP₆₀ and MP₁₈₀]. Similarly, isolated hexagons have a homogenous energy distribution over their bonds. All in all, isolated pentagons store, on average, higher energies per bond than the hexagons. When considering the whole polygons, the energies become comparable. Importantly, the stiffness of the springs affects the average bond energies (due to the interaction potential) but energy distribution remains unchanged (Fig. S6).

During the initial steps of the applied deformation, the strain reflects as an energy increase on the bonds neighboring the stress application region. In particular, the energy of the bonds aligned nearly in parallel with the direction of the exerted stress (dark red) is higher than that of the transversal bonds (dark blue). This behavior is preserved throughout the pulling simulations, and the results show that at the maximum extension, the bond energies increase symmetrically with increasing distance from the center of the MP, with the bonds closest to the pulling extremes having the highest energies [Fig. 4(c)]. Hence, the distribution of stress over the bonds is heterogeneous and varies with system size. Given a metaparticle with three links per bead, the distribution of stress is done over a maximum of 3n/2 bonds, where n is the number of beads.

To minimize the average energy per bond under high stress conditions, the metaparticles undergo internal rearrangements. The rearrangements can be spontaneous, i.e., in one step, or sequential, i.e., in two or multiple steps. The small MPs (R_g < 3σ) deform in one step along the direction of the applied stress independently of their symmetry. In this case, the excluded volume interactions between the beads prevent the system from undergoing any internal deformations. For the larger systems with D_6d symmetry, the topology largely influences which deformation mechanism the metaparticles follow. For instance, MP₁₈₀ and MP₂₄₀ follow a multi-step deformation. In the first steps of the deformation, the MP extends along the direction of the pull. Then, the MP contracts and attains a flower-like conformation determined by the arrangement of the pentagons (Figs. S8 and S9) and the direction of the pull.

Next, under continued applied stress, the complex rearranges and bends inward upon itself as response to the applied forces. This arrangement is preserved throughout the simulation or until the maximum extension limit has been reached. As eluded, the rearrangement depends on the topology of the MP and the arrangement of the pentagons relative to the direction of the pull due to the anisotropy in surface tension.

Next, we investigated the ability of the metaparticles to recover to their initial conformation upon stress release. For this, we compared the time required for the MPs to reach their maximum extensibility (evaluated at the UTS), τ^d, with the recovery time, τ^r, i.e., the time required for the MPs to return to their initial conformation after stress release (Fig. 5). A system is considered as fully recovered, when the transversal and axial diameters are within the standard deviation of their respective initial equilibrated values.

The analysis reveals that the deformation time increases with metaparticle size, except for MP₁₈₀ and MP₂₄₀, which have comparable τ^d. Furthermore, stiffer metaparticles require longer times to reach maximum extensibility (Fig. S10). In addition, the longer τ^d may be associated with the deformation mechanism, as MPs undergoing multistep deformations require longer times to reach their maximum extension.

The comparison between the deformation and recovery times shows a linear correlation for small MPs (R_g < 4σ), which require comparable times to deform and to recover independently of their topology (Fig. 5). The relaxation time of MP₄₈ is about three times longer than its deformation time due to the direction of applied stress, i.e., the forces are applied along the shortest axis due to symmetry reasons. With increasing system size, topology and metaparticle deformation affect differently the two timescales. For instance, the deformation and internal rearrangement of MP₂₄₀ reflects in the longer τ^r compared to τ^d (72.8τ and 32.1τ, respectively). The longer τ^r for MP₂₄₀ can be linked to the larger number of bonds that are involved in the relaxation process compared to other systems. Furthermore, visual inspection revealed that the complete recovery to the initial state is limited by the conformation of a pentagon, which remains in the “inward” folded configuration. This unique conformation reflects minimally on the overall conformation of the MP. Extended runs revealed that the complete recovery of this pentagon extends on 100-fold longer timescales.

Similar to the multi-step deformation, the recovery of the metaparticles can also occur following different pathways, i.e., single or multi-step relaxation. On the one hand, the relaxation of MPs that do not undergo internal rearrangements occurs in a single step, i.e., the particle contracts upon stress release. The thermal fluctuations and the FENE potential gradually drive the MP to its initial conformation. On the other hand, when MPs undergo internal rearrangements, the system requires a longer recovery time due to the multi-step recovery process. For instance, upon stress release, MP₄₈ first axially contracts, while maintaining its “folded” conformation. Then, the deformed surfaces recover one by one to their equilibrium state. The complete recovery to the initial state as a response to stress is a property of elastomeric materials, particularly bio-inspired ones.⁹⁰ 

We developed a coarse-grained model for MetaParticles, i.e., soft, resilient particles with tunable and responsive properties. This was achieved by representing a MP as a collection of beads interconnected by springs. The topological arrangement of the beads and the bonds is symmetrical and was inspired by the topologies of closed fullerenes. Hence, each bead is linked to three surrounding beads forming a pattern of hexagons and 12 pentagons, endowing the molecules with quasi-spherical shapes. This topological arrangement brings unique symmetries to each metaparticle. Our results show that the MPs are responsive to the applied stress following an elastomeric behavior driven by a combination of size, topology (symmetry) and area of applied stress. The deformation is reversible as the metaparticles return to their original shapes upon stress release. Furthermore, the induced deformation and the relaxation occur on different timescales depending on the deformation and relaxation mechanisms. Single step deformation leads to a linear correlation between the two timescales and is typical for MPs with R_g < 3 nm. Multistep deformations and relaxations require longer times due to the internal rearrangements of the beads, a feature shared by the larger metaparticles.

As discussed in Sec. I, the development of a particle with tunable and responsive properties constitutes a major step toward modeling and understanding the behavior of deformable particles in the biological environment. Our aim is to introduce a generic model, easy to fine-tune to include specificity and to be used to understand the responsiveness and adaptability of semi-flexible nanocarriers to environmental, chemical, and physical changes present in biological systems.

Different nanoparticles, such as lipid nanoparticles, protein cages, polymer-based nanocarriers, and dendrimers, possess unique structural and chemical properties that require specific adjustments to the model. Hence, depending on the system of interest, one may consider different parameterization schemes of the MPs. For nanocarriers with drug delivery purposes, we envision a model, in which the beads of the MetaParticles represent a collection of ligands on the surface of a nanocarrier meant to enhance adhesion to specific membranes and cellular uptake. In this case, one may consider first employing molecular dynamics simulations at full atomistic resolution to explore the interaction of specific ligands (e.g., peptides, proteins, and antibodies) and the cellular membrane. Furthermore, using enhanced sampling techniques, one may extract binding free energies, which then translate into effective potentials in the metaparticle model. Just like biological macromolecules, which can have different protonation states under distinct conditions, the coarse-grained beads will be able to adjust their interactions in response to the environment. For instance, they can become attractive or repulsive depending on the proximity to or embedding in the membrane, hence capturing the specific properties of the compounds at different pH. In context of protein-based nanocarriers,⁹¹ one may start from the high resolution structures of the nanocarrier of interest^92–94 to fine-tune the arrangement of the beads on the surface of the metaparticle. The flexibility of the cage may be first analyzed using elastic network models (ENMs) to determine the vibrational modes that arise from the collective motion of the proteins.⁹⁵ For small systems, short atomistic simulations may be carried out and analyzed to explore the link between the packing of the proteins and the flexibility of the complex.⁹⁶ Furthermore, a normal mode analysis using classical^97–99 or machine learning techniques^100,101 could be used to determine the dynamic properties of the nanocluster. These properties can then be incorporated in the flexibility of the coarse-grained model to potentially consider a dynamic adaptability of the bead stiffness. In addition, from quantities such as atomic mean-square displacements or B factors, the packing of the macromolecules can be accounted for in the metaparticle model by endowing the springs with distinct spring constants.

As highlighted throughout this manuscript, the generic behavior of this model extends beyond pure drug delivery understanding and aids in the design of metamaterials with bioinspired properties.¹⁰² For example, the mechanical properties of granular metamaterials made from flexible tessellations could be investigated by parameterizing the MP model based on existing experimental and simulation data.^103,104 Given the broad pool of applications, future questions will address the effects of bead arrangement, topology, defects, and activity on the mechanical properties of the metaparticles.

The supplementary material contains further analysis of the physical properties of the metaparticles, their mechanical properties upon changing the stiffness of the springs, the energy distribution per bond for all metaparticles, and the timelines of the deformation and relaxation times.

I.M.I. acknowledges support from the Sectorplan Bèta & Techniek of the Dutch Government, the Dementia Research - Synapsis Foundation Switzerland, and the Molecular Material Design Technology Hub.

The authors have no conflicts to disclose.

Massimiliano Paesani: Conceptualization (supporting); Data curation (lead); Formal analysis (equal); Investigation (lead); Methodology (equal); Software (lead); Validation (lead); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Ioana M. Ilie: Conceptualization (lead); Data curation (supporting); Formal analysis (equal); Investigation (supporting); Methodology (supporting); Project administration (lead); Resources (lead); Software (supporting); Supervision (lead); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The input files are publicly available on: https://github.com/ilieim/MetaParticle.