Hugo Le Roy
email: h.leroy@epfl.ch

M Mert Terzi

Martin Lenz
email: martin.lenz@universite-paris-saclay.fr

Collective deformation modes promote fibrous self-assembly in protein-like particles

The self-assembly of particles into organized structures is a key feature of living organisms and a major engineering challenge. While it may proceed through the binding of perfectly matched, puzzle-pieces-like particles, many other instances involve ill-fitting particles that must deform to fit together. These include some pathological proteins, which have a known propensity to form fibrous aggregates. Despite this observation, the general relationship between the individual characteristics of the particles and the overall structure of the aggregate is not understood. To elucidate it, we analytically and numerically study the self-assembly of two-dimensional, deformable ill-fitting particles. We find that moderately sticky particles tend to form equilibrium self-limited aggregates whose size is set by an elastic boundary layer associated with collective deformations that may extend over many particles. Particles with a soft internal deformation mode thus give rise to large aggregates. Besides, when the particles are incompressible, their aggregates tend to be anisotropic and fiber-like. Our results are preserved in a more complex particle model with randomly chosen elastic properties. This indicates that generic protein characteristics such as allostery and incompressibility could favor the formation of fibers in protein aggregation, and suggests design principles for artificial self-assembling structures.

Functional structures in living cells are often selfassembled from several copies of a single protein, from microtubules and clathrin cages to viral capsids in the shape of cylinders or spheres [START_REF] Amos | Arrangement of subunits in flagellar microtubules[END_REF][START_REF] Pearse | Clathrin: a unique protein associated with intracellular transfer of membrane by coated vesicles[END_REF][START_REF] Caspar | Physical principles in the construction of regular viruses[END_REF]. The radius of such assemblies is dictated by the curvatures of the individual particles that precisely fit together to form them. Similarly, artificial self-assembly often relies on fitting welladjusted particles together to build structures with a controlled size [START_REF] Wagenbauer | Gigadaltonscale shape-programmable DNA assemblies[END_REF][START_REF] Sigl | Programmable icosahedral shell system for virus trapping[END_REF][START_REF] Hayakawa | Geometrically programmed self-limited assembly of tubules using DNA origami colloids[END_REF][START_REF] Sacanna | Lock and key colloids[END_REF].

In other instances however, the shapes of the individual particles are ill-fitting and do not obviously dictate the structure of the aggregate. This is the case in the pathological aggregation of normally soluble proteins, i.e., of proteins not evolutionarily optimized to selfassemble into a well-defined structure [START_REF] Eaton | Sickle cell hemoglobin polymerization[END_REF][START_REF] Nelson | Recent atomic models of amyloid fibril structure[END_REF][START_REF] Ranson | Insights into the architecture of the Ure2p yeast protein assemblies from helical twisted fibrils[END_REF][START_REF] Elam | Amyloid-like filaments and water-filled nanotubes formed by SOD1 mutant proteins linked to familial ALS[END_REF]. Despite the diversity of the shapes and interactions involved, the aggregation of these ill-fitting proteins produces fibrous structures with remarkable consistency. These fibers display varied widths and internal structures [START_REF] Bousset | Structural and functional characterization of two alpha-synuclein strains[END_REF][START_REF] Fitzpatrick | Cryo-EM structures of tau filaments from Alzheimer's disease[END_REF][START_REF] Alam | αsynuclein oligomers and fibrils: a spectrum of species, a spectrum of toxicities[END_REF], and the proteins within are often significantly deformed in ways that depend on the assembly protocol [START_REF] Eichner | A diversity of assembly mechanisms of a generic amyloid fold[END_REF]. Deformations are common in proteins, and many display physiologically relevant deformation modes that facilitate self-assembly [START_REF] Perotti | Elasticity theory of the maturation of viral capsids[END_REF][START_REF] Bruinsma | Physics of viral shells[END_REF], perform a motor function [START_REF] Chappie | A pseudoatomic model of the dynamin polymer identifies a hydrolysis-dependent powerstroke[END_REF], participate in their biochemical activity [START_REF] Mapa | The Conformational Dynamics of the Mitochondrial Hsp70 Chaperone[END_REF], or serve to mechanically transmit a signal, a function known as allostery [START_REF] Wodak | Allostery in Its Many Disguises: From Theory to Applications[END_REF][START_REF] Hayer-Hartl | The GroEL-GroES Chaperonin Machine: A Nano-Cage for Protein Folding[END_REF]. Nevertheless, the generic implications of the deformability of proteins on their ill-fitting aggregation is not understood.

Beyond proteins, particle deformations have long been suggested as a mechanism to regulate aggregate size in self-assembly [START_REF] Penrose | Self-reproducing machines[END_REF]. In this picture, ill-fitting particles are forced to deform as they tightly bind to one another. As more and more particles are added to the aggregate, the distortions build up until they become so severe as to prevent any further assembly [START_REF] Hagan | Equilibrium mechanisms of self-limiting assembly[END_REF][START_REF] Spivack | Stress accumulation versus shape flattening in frustrated, warpedjigsaw particle assemblies[END_REF]. The accumulation of stresses resulting from such distortions may govern the structure of DNA origami assemblies [START_REF] Berengut | Self-limiting polymerization of DNA origami subunits with strain accumulation[END_REF][START_REF] Videbaek | Tiling a tubule: how increasing complexity improves the yield of self-limited assembly[END_REF] and prevent the indefinite bundling of preexisting protein fibers [START_REF] Claessens | Helical twist controls the thickness of Factin bundles[END_REF][START_REF] Adamcik | Understanding amyloid aggregation by statistical analysis of atomic force microscopy images[END_REF][START_REF] Brown | An equilibrium double-twist model for the radial structure of collagen fibrils[END_REF]. Beyond merely fixing the overall size of an aggregate, the accumulation of deformations can moreover dramatically alter its shape. This has been proposed to drive a transition from cylindrical to tapelike fiber bundles [START_REF] Hall | Morphology selection via geometric frustration in chiral filament bundles[END_REF][START_REF] Hall | How geometric frustration shapes twisted fibres, inside and out: Competing morphologies of chiral filament assembly[END_REF]. Finally, it can also drive sticky, deformable particles to form anisotropic aggregates that grow into infinite one-dimensional structures reminiscent of pathological protein fibers [START_REF] Lenz | Geometrical frustration yields fibre formation in self-assembly[END_REF]. The underlying mechanism and the nature of the particle properties that determine the dimensionality of the final aggregate however remain elusive.

In this paper, we provide a detailed analytical understanding of the emergence of self-limited and fibrous aggregates in a two-dimensional deformable particle system. We first introduce a minimal model based on highly symmetrical particles. It gives rise to an emergent elastic boundary layer length ℓ, allowing us to map it onto a continuum description in the limit of large ℓ. We use this description to compare the energies of several candidate structures and establish an aggregation phase diagram, which we then validate using numerical simulations. Finally, we introduce a much broader class of elastic particles, and demonstrate that the results derived in the idealized model still apply there, including in cases where the values of ℓ are moderate.

FIG. 1. The assembly of ill-fitting particles results in collective deformations, shown here in a minimal 1D model. (a) left: Individual particle at rest, right: colored springs schematizing the elasticity of the particles. In this 1D model, particles are allowed to aggregate only along the horizontal direction. Their corners then become colocalized (black arrows), which requires deforming at least some of the springs. (b) Schematic of a one-dimensional particle aggregate showing the state of the springs therein. While the yellow and red springs are able to assume different lengths in the vicinity of the edge of the aggregate, they are forced to have the same lengths in the bulk. The resulting energetic penalty hampers the formation of space-filling, bulky aggregates. In two or three dimensions, a similar penalty may result in the formation of fibrous aggregates.

ELASTIC AGGREGATION MODEL

To understand the interplay between particle deformations and aggregate structure, we first discuss a minimal one-dimensional example where ill-fitting particles deform upon aggregation. This leads to a deformation gradient from moderately deformed particles at the edge of the aggregate to highly deformed ones in its bulk. We then introduce a two-dimensional model that allows for a much wider diversity of aggregate structures, including fibers and planes. This model is analytically intractable in its general form, which leads us to design a continuum limit to enable further analysis.

One-dimensional toy model

We consider a collection of identical isosceles trapezoids [Fig. 1(a)]. Each such particle can aggregate with its left and right neighbors by fusing its vertical sides with theirs. In the special case where the trapezoids are welladjusted, i.e., if they are rectangles, such binding does not require any deformation. Conversely, particles whose top and bottom faces have different lengths are ill-fitting and must deform to bind. We model the energetic cost of this deformation using four springs: two representing the top and bottom faces of the particles (yellow and red) with rest lengths 1 ± ϵ and spring constant k, and two connecting springs with spring constants k c /2 and rest lengths ϵ that tend to center the top and bottom faces (blue). Summing the contributions of these four springs, we write the deformation energy of the central particle of Fig. 1(b) as

e (i) d = k 2 [(x ↑ i+1 -x ↑ i ) -(1 + ϵ)] 2 + k 2 [(x ↓ i+1 -x ↓ i ) -(1 -ϵ)] 2 + k c 4 (x ↑ i+1 -x ↓ i+1 -ϵ) 2 + k c 4 (x ↓ i -x ↑ i -ϵ) 2 , (1)
where the {x ↑ i }, {x ↓ i } denote the coordinates of particle corners. This model does not involve any explicit prestresses; including some would not make any difference within the linear response regimes studied here and in the rest of this work.

Defining the shift between an upper and lower corner as

δ i = x ↑ i -x ↓ i , force balance dictates that inside the aggregate k(δ i+1 -2δ i + δ i-1 ) = k c δ i ⇒ δ i ∝ ϵ sinh(i/ℓ), (2)
where we define i = 0 as the center of the aggregate and where

ℓ = 1/ ln 1 + k c /k + 2k c /k + (k c /k) 2 ∼ kc≪k k/2k c . (3)
The full prefactor of the last expression of ( 2) is fixed through the force balance condition at the aggregate's left and right edges. In a large aggregate, it results in an exponential decay δ i ∝ ϵ exp(-|ii edge |/ℓ) close to these edges. The initially trapezoidal particles at the center of the aggregate are thus forced into a rectangular shape (δ i = 0), in contrast with the particles that reside within an edge-associated elastic boundary layer of size ℓ.

In the limit k c /k → 0, the boundary layer size diverges. This regime is characterized by very rigid yellow and red springs, implying that the yellow and red springs close to the edge of the aggregate are almost at their equilibrium lengths. Going deeper into the aggregate, each blue spring exerts a small compressive (tensile) force on the yellow (red) chain. These forces add up over long distances, implying a progressive change of the yellow and red strain over a length scale much larger than the particle size.

Two-dimensional particle model

One-dimensional aggregates are very simple geometrically, and are entirely characterized by the number of particles that they contain. Aggregation in higher dimension allows for a much broader variety of aggregate structures. To study the emergence of complex shapes there, we introduce a two-dimensional model based on hexagonal particles. Well-adjusted particles are represented by regular hexagons. Ill-fitting hexagons, by contrast, have alternating protruding (yellow) and withdrawn (red) corners [Fig. 2(a)]. The corners belonging to each category form a yellow and a red equilateral triangle whose sides are springs with rest lengths 1 ± ϵ and a spring constant k. The (blue) sides of the hexagon itself play the role of connecting springs with spring constant k c , and are at rest when the yellow and red equilateral triangles are themselves at rest and centered. These hexagonal particles are three-fold symmetric, which rules out an intrinsic, particle-level preference for forming one-dimensional fibers. Any fiber formed from their aggregation will thus be an emergent symmetry-broken structure. The response of these particles to shear and uniform compression cannot be independently varied while holding k c /k (and therefore the boundary layer size) constant. To enable particles that range from fully compressible to incompressible, we thus additionally endow both yellow and red triangles with an areal rigidity. We implement it through an energy e

↑/↓ area = k area (A ↑/↓ -A ↑/↓ 0 ) 2 /2A ↑/↓ 0 ,
where ↑ and ↓ respectively refer to the yellow and red triangle and A (A 0 ) are the associated triangle areas (rest areas).

The quadratic spring and areal energies introduced above all vanish in the particle's resting state. Any deformation away from this state implies an energetic cost, and such deformations are required to accommodate particle binding. In our model, two particles can bind along a blue side by merging one yellow and one red corner each. The merging of corners with different colors is not allowed. Each pair of bound sides is rewarded by an energy -g regardless of the particles' state of deformation, which defines a zero-range interaction between particles. These rules favor the assembly of hexagons into a triangular "particle lattice" [Fig. 2(b)] where all pairs of neighboring particles are bound, which we consider throughout. The aggregate topology, i.e., the specification of which particles bind to which others through which sides, can thus be entirely described by considering a triangular lattice and specifying a list of the lattice sites that are occupied by a particle. In the following we use the symbol T to refer to this topology. Since the binding energy is fully determined by the number of bound particle sides, it only depends on T .

Continuum formalism

Finding the most favorable aggregate in our 2D particle model requires two steps: to compute the optimal deformation energy for each fixed topology T , and then to determine which topology has the lowest optimal energy. Here we introduce a continuum approximation that renders the first step analytically tractable in several important cases. This approximation is formally valid in the limit where the 2D counterpart of the boundary layer size ℓ is much larger than the particle size.

To define our continuum limit, we note that in large aggregates where all sites of the particle lattice discussed above are occupied (i.e., without holes), the yellow and red springs arrange into triangular spring lattices. In the regime k c ≪ k where connecting (blue) springs are much softer than triangle (yellow and red) springs, the strain within the yellow and red triangular spring lattices varies slowly over space. This is similar to the behavior of our one-dimensional model. As a result, we can assimilate each of these triangular spring lattices to a continuum sheet, giving rise to a continuum elastic energy

E d =   λ 2 ∂ α u ↑ α -2ϵ 2 + µ ∂ α u ↑ β + ∂ β u ↑ α 2 -ϵδ αβ 2   dA +   λ 2 ∂ α u ↓ α + 2ϵ 2 + µ ∂ α u ↓ β + ∂ β u ↓ α 2 + ϵδ αβ 2   dA + κ c 2 (u ↑ α -u ↓ α ) 2 dA, (4)
where the superscripts ↑ and ↓ refer to the yellow and red sheets respectively, and where the summation over repeated indices is implied while δ denotes the Kronecker delta. The displacement fields u ↑/↓ (r) of either sheet are computed with respect to the infinite-aggregate, bulk state where all hexagons are regular, a state akin to a row of length-one rectangles in the one-dimensional model. Neither elastic sheet is at rest in this reference state, and r is the position vector in this state. The displacement gradient ∂ α u ↑/↓ α thus plays the same role as the finite difference (x [START_REF] Amos | Arrangement of subunits in flagellar microtubules[END_REF]. The first integral of ( 4) is a two-dimensional generalization of the first term of (1), and gives the elastic energy of an isotropic elastic sheet with Lamé coefficients λ and µ whose resting state is characterized by an isotropic strain ∂ α u ↑ β = ϵδ αβ . As a reminder, µ is the usual shear modulus and the sheet's Young modulus and Poisson ratio are Y = 4µ(λ + µ)/(λ + 2µ) and ν = λ/(λ + 2µ). The integration element dA runs over the reference area A. Finally, the last integral captures the energy of the connecting springs. It penalize any shift between the centers of the yellow and red triangles of a given particle, and therefore any difference in the displacements of the two sheets.

↑/↓ i+1 -x ↑/↓ i -1) of
The harmonic form of ( 4) is valid for small deformations, which can be obtained for any aggregate topology considered here given a small enough ϵ. Within this limit, the particles' binding energy can be described through

E b = γL, (5)
where γ is a line tension and L is the total length of aggregate edge in the reference state, including any internal holes. The parameters of the continuum model are mapped onto those of the 2D particles in the Supporting Information, yielding

µ = √ 3 4 k (6a) ν = √ 3k + 2k area 3 √ 3k + 2k area (6b) κ c = 2 √ 3k c (6c) γ = 4 √ 3 g (6d)
The binding energy of an aggregate depends only on topology T , which we denote as E b (T ). By contrast, the deformation energy E d (T , {u ↑/↓ (r)}) depends both on the topology and on the displacement fields u ↑/↓ (r). As described in the beginning of this section, the optimal energy associated with a given aggregate topology is thus

E(T ) = E b (T ) + min {u ↑/↓ } E d (T , {u ↑/↓ }). (7)
Once this minimization is performed, finding the most favorable aggregate structure requires finding the topology T that minimizes E(T ).

AGGREGATION PHASE DIAGRAM

To establish an aggregation phase diagram, we consider a system with a fixed but large number of particles and ask which binding topology minimizes the total energy of the system. We first use our continuum formalism to compute the energies of an infinite bulk, an elongated fiber and a disk-like aggregate, thus offering a first comparison of the stability of two-, one-and zero-dimensional structures. We then numerically compare the energies of a wider range of putative aggregate structures in our discrete particle model. Finally, we confirm the converging results of these two approaches using numerical simulations devoid of a priori constraints on the aggregate topology.

Continuum phase diagram

In 2D space-filling, infinite aggregate, all particles are forced into a regular hexagonal shape. The continuum energy of (4) is then minimal for u ↑ = u ↓ = 0, yielding an optimal energy per unit reference surface

e bulk = µ 1 + ν 1 -ν 4ϵ 2 . ( 8
)
Denoting the reference area per particle by a, the total energy of a set of N particles thus reads N ae bulk in the large-N limit. The binding energy is proportional to the perimeter of the aggregate (E b ∝ √ N ), and is thus negligible.

To determine whether fiber formation is favored over bulk aggregation, we minimize (4) over {u ↑/↓ (r)} for an infinite strip of width W and find (Supporting Information)

u ↑ (x) = -u ↓ (x) = ℓϵ (1 + ν) sinh(x/ℓ) cosh(W/2ℓ) x, ( 9
)
where x is the direction perpendicular to the fiber and where

ℓ = λ + 2µ 2κ c . ( 10
)
The quantity λ + 2µ is known as the P-wave modulus of the sheet, and charcterizes the cost of compressing it along one axis without allowing it to deform in the perpendicular direction. Similar to (1) and Fig. 1(b), the profile of ( 9) implies bulk-like, highly deformed particles in the center of the fiber, while close to the edges the red and yellow sheets gradually relax within a boundary layer of width ℓ. Defining the dimensionless line tension Γ = 2γ/[(1 + ν)ℓe bulk ], line tension cost associated with the fibers' edges reads E b = 2ΓN a/W and the mean energy per unit surface reads e fiber (W )

e bulk = 1 -(1 + ν) tanh(W/2ℓ) W/ℓ + (1 + ν) Γ W/ℓ . (11)
In the ϵ → 0, small-particle-mismatch limit, all deformation energies scale as ϵ 2 . Thus the parameter Γ encloses both the γ and the ϵ dependence of all self-assembly outcomes studied in this paper. When Γ < 1, e fiber (W ) displays a minimum at a finite fiber width W * . This optimal width diverges in the limit Γ → 1, and the corresponding fiber is always more stable than the bulk [Fig. 3(a)]. To understand this stability, consider a semi-infinite aggregate that fills half of the plane. While its energy per unit surface far from its edge is equal to e bulk , the presence of the edge brings about two energetic contributions. The first is a bare line tension cost γ per unit edge length. The second is the deformation energy gain in the boundary layer, which is of the order of e bulk per unit area. Since the width of the layer is ℓ, the resulting energy gain per unit edge length is of the order of ℓe bulk . For γ ≲ ℓe bulk , forming a new edge thus results in a net energy gain.

At the scaling analysis level, this is equivalent to Γ < 1. This argument implies that infinite bulks can lower their energy by breaking up into fibers in this regime. However, if these fibers are made so narrow that their widths become of order ℓ or smaller, the boundary layers associated with their two edges start to overlap. Such narrow fibers can only claim a fraction of the deformation energy reduction described above. As a result, very narrow fibers are penalized. This implies the existence of an optimal width W * that is of order ℓ when Γ is of order one but smaller than one.

In an aggregate whose resting shape is a disk of radius R, the displacement field is given by (Supporting Information)

u ↑ (r) = ℓϵ (1 + ν)I 1 (r/ℓ) I 0 (R/ℓ) + I 2 (R/ℓ) + ν[I 0 (R/ℓ) -I 2 (R/ℓ)] r (12
) and u ↓ (r) = -u ↑ (r). Here r denotes the radial coordinate and the I α s are the modified Bessel functions of the first kind. Just like Eqs. ( 2) and ( 9), ( 12) predicts an exponential decay of the boundary displacement at the edge of large aggregates. Again, the deformation cost is smaller at the aggregate edge than in its center, yielding an average energy per unit area

e disk (R) e bulk = 1 - I 0 (R/ℓ) -I 2 (R/ℓ) I 0 (R/ℓ) + I 2 (R/ℓ) + ν[I 0 (R/ℓ) -I 2 (R/ℓ)] + (1 + ν) Γ R/ℓ . ( 13
)
This expression displays an optimal finite aggregate size R * at low values of Γ. As in the fiber case, this optimal disk is more stable than the bulk up to values of Γ of order one, although the exact criterion differs due to the curved geometry of the interface [Fig. 3(a)].

Our phase diagram indicates that fibers are more stable than disks for large Poisson ratios, i.e., they are favored in the aggregation of incompressible particles (characterized by ν = 1 in 2D). To understand this, we compare a vertical fiber and a disk at ν = 1. Symmetry forbids vertical (orthoradial) displacements in the fiber (disk), allowing only horizontal (radial) displacements. However, no such displacement is possible without violating the incompressibility condition of the yellow or red sheet. As a result, the sheets must remain in their resting states, implying displacements u

↑ x = -u ↓ x = 2ϵx (u ↑ r = -u ↓ r = 2ϵr
) to lowest order in ϵ with respect to the fictitious bulk reference state. All the deformation cost thus comes from the connecting springs. The resulting connecting spring energy per unit area reads 8κϵ 2 x 2 (8κϵ 2 r 2 ), proportional to the square of the distance from the center of the aggregate. This is where the difference between fibers and disks manifests itself. The edge of a fiber is just as long as its centerline, while the center of a disk is much smaller than its perimeter. As a result, a smaller proportion of connecting springs are highly extended in the fiber than in the disk, making the former energetically cheaper. Conversely, in the limit of low Poisson ratio ν → 0, the gain per unit length of the straight, fiber-like and of the curved, disk-like boundary layers become identical. Forming such boundary layers is favorable overall for Γ < 1. Since the disk tends to have more boundary layer per unit area than the fiber, it is more favorable in this limit.

Comparison of pre-made discrete aggregate structures

The phase diagram of Fig. 3(a) focuses on continuum sheets, leaving open the question of whether the formation of holes within the aggregate or regimes where the particle size is comparable to the boundary layer thickness could result in different aggregation behaviors. To assess its robustness to these effects, we numerically implement our discrete particle model in a computer. We consider several aggregates with predetermined topologies including periodic bulks with and without holes, fibers of various widths as well as hexagonal aggregates approximating disks of different radii [Fig. 4(a)]. Isolated particles are taken into account as disks with one particle. We use a conjugate gradient algorithm to minimize the energy of each aggregate over the position of the all particle corners, analogous to the minimization over the deformation field in [START_REF] Sacanna | Lock and key colloids[END_REF]. We build three different fiber structures by cutting the bulk along distinct directions. One of the fibers is left-right asymmetric and spontaneously curves, although that curvature vanishes in the large-W limit. In practice we obtain the energy of infinitely long fibers by extrapolating from long ones with increasing lengths.

To compute the aggregation phase diagram, we minimize the deformation energy of aggregates with width (radii) ranging from 1 to 35 [START_REF] Berengut | Self-limiting polymerization of DNA origami subunits with strain accumulation[END_REF] for ℓ = 5 and for values of ν ranging from 0 to 1. For each value of Γ in Fig. 3(b) we select the aggregate with the lowest total energy. Bulks with holes, which in our model have zero deformation energy, dominate the assembly only at very small surface tensions. The rest of the phase diagram is essentially identical to the continuum one, except for an expansion of the fiber region against both bulks and disks. This may be due to the increased stability of curved fibers, which are always more stable than their straight counterparts.

To assess the consistency of the morphology of finiteℓ aggregates with the continuum (ℓ → ∞) expectation of Fig. 3(a), we numerically determine the most favorable fiber width and disk radius for a range of ℓ and Γ in Fig. 4(b). While the two approaches are guaranteed to agree only in the large-ℓ limit, in practice the continuum approximation yields very accurate predictions all FIG. 4. Comparison of the energies of a collection of aggregates with pre-determined topologies. (a) List of the aggregate topologies included in the trial: bulk with holes, three types of fibers obtained by piling the particles in three different ways, hexagonal aggregates (including single particles) and bulk. The white lines and shading in the bulk show the three different types of cuts used to produce the fibers shown on the left. In practice, the isolated particles and uncurved fibers are never the most stable structures. We compare a broad range of fiber widths and hexagonal aggregates, despite representing only a few here. (b) left: When measured in units of number of particles, the optimal fiber width W * (computed here for ν = 0.95) and disk radius R * (ν = 0.2), strongly depend on the particles' elastic properties through ℓ. Right: Rescaling both lengths by ℓ however leads to an excellent collapse with the analytical prediction (dashed red line), even for small ℓ.

the way down to values of ℓ equal to the particle size. This indicates that the boundary layer physics revealed by our continuum model remains an excellent qualitative and quantitative description of the aggregation process even when the stiffness of the connecting springs k c is comparable to that of the others (k c ≲ k ⇔ ℓ ≳ 1). The length ℓ thus provides a robust tool to predict the typical number of particles in the cross-section of a fiber or a disk.

Monte-Carlo validation of the phase diagram

As a final validation of our phase diagram, we remove any restriction on the aggregate's structure and use a Monte-Carlo algorithm to evolve its topology. We simulate a triangular lattice where each site can be empty or occupied by a particle. We start with randomly placed particles, and attempt Monte-Carlo moves where a ran-domly chosen particle is moved to a randomly chosen empty site. We compute the optimal energy of the resulting new topology using our conjugate gradient method. The move is accepted according to a Metropolis criterion with temperature T . Since the deformation energy is optimized before the application of the Metropolis criterion, this temperature applies only to the system's topological degrees of freedom. To look for an approximation of the system's topological ground state, we perform a simulated annealing procedure whereby T is slowly lowered from a large value to zero over the course of the simulation.

The computational cost and limited particle number in our simulation make it difficult to construct a full Monte-Carlo phase diagram. Instead, we simulate select parameter regimes to validate our main findings. Specifically, we perform three line scans at ℓ = 5 whose final aggregation states are shown in the small panels of Fig. 3(b). First, a range of increasing Γ at ν = 0.9 shows the dominance of fibers at large Poisson ratios, and the overall tendency of the fiber widths to increase with increasing line tension. A second horizontal scan at ν = 0.4 shows disk-like aggregates whose radii increase with increasing Γ. Finally, a vertical scan at fixed Γ = 0.5 shows a transition between disks and elongated aggregates at the predicted value of ν. The shapes of the aggregates resulting from the simulations are somewhat variable, as can be assessed from the Γ = 0.5, ν = 0.9 condition which belongs to two different line scans and for which the outcomes of two independent simulations is shown in Fig. 3(b). Our results show that despite this limitation and the relatively small number of particles in our system, our phase diagram accurately predicts the overall outcome of an unconstrained assembly.

EXTENSION TO RANDOM PARTICLES

The discrete model of Fig. 2 has been specifically designed with the continuum approximation in mind, making it unclear whether our results apply to more generic particle types. To test whether this is the case, we define a much broader class of hexagonal particles with the same shape at rest and the same binding rules but with a considerably more generic deformation energy. We parametrize the shape of each particle by nine distances d 1 , d 2 , . . . , d 9 defined in Fig. 5(a). The vector d of these lengths completely characterizes a particle's shape, and we denote its value at rest by d 0 . The deformation energy of one particle is an arbitrary quadratic form of the deviation from that state:

e d = 1 2 (d -d 0 ) T • M • (d -d 0 ). ( 14
)
In the following we draw the matrix M from a random distribution that ensures that it is positive semidefinite and three-fold symmetric (Supporting Information). Such particles allow not only for differences in the elastic constants of the colored springs of Fig. 2(a), but also for new couplings between them. For instance, one of these new couplings dictates that compressing a yellow spring makes the red spring to its right shrink or extend. It also opens regimes where the coupling springs are not much softer than the others. We first determine whether the length ℓ in a random particle is determined by the same ratio of uniaxialcompression to triangle-shifting moduli as in [START_REF] Ranson | Insights into the architecture of the Ure2p yeast protein assemblies from helical twisted fibrils[END_REF]. To access a wide range of ℓ values, we randomly draw a large number (3 × 10 7 ) of instances of matrix M. For each instance, we operationally define ℓ as the thickness of the elastic boundary layer in a simulation of a semi-infinite aggregate [Fig. 5(b)]. We then relate this measured value to the properties of a single particle, by numerically measuring random-particle proxies of the two moduli λ + 2µ and κ c and using [START_REF] Ranson | Insights into the architecture of the Ure2p yeast protein assemblies from helical twisted fibrils[END_REF] to compute a predicted boundary layer thickness ℓ pred (Supporting Information). As shown in Fig. 5(c), we find a good correlation between ℓ and ℓ pred for a randomly selected subsample of 10 5 matrices. This indicates that frustration at the edge of an aggregate relaxes in similar ways in simple and random particles, with the two subtriangles of the particles progressively shifting relative to each other. Similar to the one-dimensional model of Fig. 1, this shift causes a restoring force (denoted k c δ i in the 1D model) which is balanced by stiffness of the particles' other degrees of freedom (associated with the stiffness k in the 1D model). Equation [START_REF] Ranson | Insights into the architecture of the Ure2p yeast protein assemblies from helical twisted fibrils[END_REF], which we use to compute ℓ pred , uses the P-wave modulus λ+2µ as a proxy for this stiffness. While this value is as an upper bound for the cost of deforming the particles, in practice our random particles may deform in more complex, less costly ways. Such deformations would result in a lower effective stiffness, and could explain why ℓ pred tends to overestimate the actual value of ℓ.

To determine whether the boundary layer thickness ℓ controls the aggregate size in the same way as in our simple model, we study a randomly selected subsample of 10 3 matrices M from our total sample. We use the same discrete aggregate procedures as in Fig. 4(b) to determine which fiber/disk width is the most favorable for Γ ranging from 0 to Γ max , where we define Γ max as the critical line tension where the bulk becomes the most favorable structure. Both fibers and disks tend to become larger with increasing Γ despite a large dispersion of the data similar to that observed in Fig. 4(b) (Supporting Information). We however show in Fig. 5(d) that just as in the case of simple particles, this dispersion is largely abolished by rescaling the aggregate size by the boundary layer size ℓ. This demonstrates that the boundary layer still controls the physics of the assembly in the random particle case. Moreover, the resulting radii and width distribution are clearly centered around the values of W * and R * predicted by the continuum limit.

We finally assess the applicability of the phase diagrams of Fig. 3 to random particles by correlating the most favorable aggregate type with a suitable measure of the "individual sheet Poisson ratio" of the analytical theory. We define this measure by decoupling the "yellow" and "red" subtriangles of our new particles from each other (i.e., set

M ij = 0 for all (i, j) / ∈ [1, 3] 2 ∪ [4, 6]
2 ), composing a lattice out of each and numerically computing their separate Poisson ratios. We use the average of these two values as our ν. We then segregate our 3 × 10 7 instances of the matrix M into three groups with ℓ < 1, 1 < ℓ < 1.5 and 1.5 < ℓ. We randomly select a few thousand particles within each group to obtain three near-uniform distribution of Poisson ratios. Fig. 6 shows the type and width of the aggregates obtained in each case. Fibers form at large ν and Γ in all three groups, thus demonstrating the robust influence of these parameters. Disks also form in the expected parameter regimes, although they tend to be replaced by bulks with holes for the smallest values of ℓ. Finally, fibers are even more predominant here than in our initial model, illustrating the broad relevance of our fiber formation mechanism upon ill-fitting self-assembly.

DISCUSSION

Protein aggregation in disease involves the selfassembly of deformable, ill-fitting sticky particles. Compact aggregates of such particles display a deformation gradient between a relatively unconstrained edge and a strongly deformed, frustrated core. This deformed core is energetically costly, while the aggregate's surface tension implies that its edge is also costly. As a result, the cheapest part of the aggregate lies in between the core and the edge, in the shallow bulk region. The minimization of the aggregate's energy thus requires its structure to comprise as much of this shallow bulk as possible while keeping both core and surface small. Here we show that for a rather general model of deformable particles, this results in aggregate size limitation and emergent anisotropy. This geometrically nontrivial optimization is a priori strongly dependent on the elastic properties of the particles. We nevertheless identify two surprisingly simple particle-level predictors of its outcome, namely an elastic screening length ℓ and particle incompressibility.

In our model, the shallow bulk manifests as an elastic boundary layer with size ℓ. This is reminiscent of the size limitation mechanisms at work in ribbons self-assembled out of planar materials with an intrinsic negative Gaussian curvature [START_REF] Aggeli | Hierarchical self-assembly of chiral rod-like molecules as a model for peptide beta-sheet tapes, ribbons, fibrils, and fibers[END_REF] or specially designed warped jigsaw puzzle particles [START_REF] Spivack | Stress accumulation versus shape flattening in frustrated, warpedjigsaw particle assemblies[END_REF]. Particles with a frustrated continuum spin-like degree of freedom and coupled incommensurate lattices in the absence of phase slips also give rise to a boundary layer [START_REF] Meiri | Cumulative geometric frustration and superextensive energy scaling in a nonlinear classical xy-spin model[END_REF][START_REF] Hackney | Dispersed, condensed and self-limiting states of geometrically frustrated assembly[END_REF][START_REF] Bak | Commensurate phases, incommensurate phases and the devil's staircase[END_REF]. In contrast with these specific particle geometries however, here this behavior emerges in a wide variety of particles with randomly chosen elastic properties. This suggests that the mechanism may also apply in packings of complex ill-fitting proteins. Our observation of a connection between soft deformation modes at the particle level, large values of ℓ and consequently large self-limited aggregates thus suggests an analogy with allosteric proteins, which mechanically transmit a signal by undergoing a concerted conformational change along a soft deformation mode [START_REF] Thirumalai | Symmetry, Rigidity, and Allosteric Signaling: From Monomeric Proteins to Molecular Machines[END_REF][START_REF] Bray | Conformational spread: the propagation of allosteric states in large multiprotein complexes[END_REF][START_REF] Bray | The propagation of allosteric states in large multiprotein complexes[END_REF].

While extended boundary layers appear to require some particle deformation modes to be much softer than others [Fig. 5(c)], not all particle-level soft modes result in one (Fig. S5). Indeed, once embedded in an aggregate, an initially soft deformation mode may couple to and be stiffened by the presence of neighboring particles. While our six-vertices, two-dimensional particles comprise only relatively simple soft modes compatible with collective deformations [namely the triangle-shifts illustrated in Fig. S4(a-b)], more complex objects are likely to allow many more. For instance, generalizations of our model in three dimensions, where the qualitative physics FIG. 6. Random particle aggregation diagram showing good agreement with the analytical results of Fig. 3. Color coding is as in Fig. 3(b). Each horizontal line in the diagrams corresponds to an instance of the elasticity matrix M. Although the boundaries between the different regions of the diagram fluctuate due to the random origin of M, together they outline very consistent regions where bulk with holes, disks, fibers and bulks dominate. The diagrams with larger values of ℓ cover a more restricted range of ν due to the relative scarcity of particles with both large ℓ and small ν in matrices M produced by our random generation procedure. See Supporting Information for details.

highlighted here still applies, could additionally involve frustration and soft modes associated with chiral particle twisting. Further investigations are required to identify the geometrical requirements for a soft mode to be compatible with collective deformations, which in turn allows the build-up of a thick boundary layer. These requirements will shed light on the specific elastic parameters that most influence the particles' aggregation behavior, which could include the formation of clusters, fibers or sheets in three dimensions. In the language of proteins, this would allow us to assess the likelihood for a given allosteric mode to dictate the size of a frustrated protein aggregate. This would constitute a very general tool for predicting the structure of an aggregate from the individual properties of its proteins.

The three-fold symmetry of our model particles implies that their propensity to form fibers is an emergent property as opposed to an intrinsic preference for uniaxial aggregation. This breaking of symmetry is reminiscent of the strain-induced, elongated structures formed during frustrated epitaxial growth [START_REF] Tersoff | Shape transition in growth of strained islands: Spontaneous formation of quantum wires[END_REF][START_REF] Spencer | Morphological instability in epitaxially strained dislocation-free solid films[END_REF][START_REF] Eaglesham | Dislocation-free Stranski-Krastanow growth of Ge on Si(100)[END_REF]. Both our most simple model and our generic, random-elasticity particles indicate that fiber formation is most advantageous in incompressible particles. This behavior is not specific to hexagonal particles, and we show in the Supporting Information that an alternative model based on triangular objects displays an essentially identical phase diagram. These results suggests that our mechanism is more relevant for dense globular proteins, which tend to be largely incompressible, than for proteins with fluffy intrinsically disordered domains. Globular proteins are indeed known to form fibers in, e.g., sickle cell anemia and amyotrophic lateral sclerosis [START_REF] Eaton | Sickle cell hemoglobin polymerization[END_REF][START_REF] Ranson | Insights into the architecture of the Ure2p yeast protein assemblies from helical twisted fibrils[END_REF][START_REF] Elam | Amyloid-like filaments and water-filled nanotubes formed by SOD1 mutant proteins linked to familial ALS[END_REF]. The binding free energies involved in these assemblies are typically at least one order of magnitude larger than the thermal energy k B T , justifying the analogy with our zero-temperature model. The particle deformations there moreover remain modest, consistent with the analytically tractable, smallfrustration (ϵ ≪ 1) regime studied here. In more asymmetric particles than those studied here, the mechanism outlined in our work could be complemented by fiberformation mechanisms based, e.g., on the presence of two specific binding sites on either side of the particle. Kinetic effects such as diffusion-limited aggregation, which hampers the formation of bulky aggregates, may also favor fibers.

Beyond proteins, the principles outlined here could be harnessed to control the assembly of artificial nanoobjects. In DNA origami, soft deformation modes can be engineered to control aggregate size in a simple onedimensional chain [START_REF] Berengut | Self-limiting polymerization of DNA origami subunits with strain accumulation[END_REF]. Two-or three-dimensional extensions of such designs should be prone to frustrationinduced fibrillation. Fibrous morphologies also emerge in DNA origami systems into which this feature is not intentionally designed [START_REF] Tikhomirov | Programmable disorder in random DNA tilings[END_REF]. While systems involving rigid nanoparticles are less straightforwardly mapped onto an elastic continuum than DNA origami, fibrous morphologies have been observed in packings of tetrahedral particles [START_REF] Yan | Self-assembly of chiral nanoparticles into semiconductor helices with tunable near-infrared optical activity[END_REF] and successfully rationalized with an elastic model [START_REF] Serafin | Frustrated self-assembly of non-euclidean crystals of nanoparticles[END_REF] based on a elastic frustration originating from a metric incompatibility. Our study does not rely on this very strong, somewhat specialized type of frustration and thus demonstrates that it is not required for frustrationinduced fiber formation. Rigid colloids with short-range attractive and long-range repulsive interactions also display frustration-induced fibrous structures [START_REF] Sciortino | One-Dimensional Cluster Growth and Branching Gels in Colloidal Systems with Short-Range Depletion Attraction and Screened Electrostatic Repulsion[END_REF]. In this case, the importance of the distance dependence of the particle interaction profile falls outside of the scope of our small-frustration (ϵ ≪ 1) formalism. More generally, frustration build-up in the presence of nonlinear elasticity and strain-induced particle unbinding can lead to the emergence of new aggregate patterns that have only begun to be explored [START_REF] Hall | Building blocks of non-Euclidean ribbons: size-controlled selfassembly via discrete frustrated particles[END_REF] and could play a crucial role in the physical implementation of the principles described here.

where additional factor of 1/2 comes from the fact that only half of the triangles are endowed with an areal stiffness in our model. Adding the two contributions of Eqs. ( S3) and (S5), we find a total sheet energy density of the form (S1) with

µ = √ 3 4 k, ν = √ 3k + 2k area 3 √ 3k + 2k area . (S6)
The second contribution to the total energy stems from the coupling springs. In the discrete particle model, the total coupling energy reads

F c = k c 2 coupling springs (r -r 0 ) 2 , ( S7
)
where r and r 0 are the deformed and rest length of the coupling springs, respectively. To lowest order in displacement, the change of length of a spring whose direction is given by the unit vector ŝ reads

r -r 0 ∼ u ↑ -u ↓ • ŝ = |u ↑ -u ↓ | cos(θ u -θ s ), ( S8
)
where θ u and θ s are the angles that the vectors u ↑u ↓ and ŝ respectively make with the horizontal axis. As required for the continuum limit approach, we assume that the displacement fields u ↑ , u ↓ are homogeneous. This yields a total energy per spring

F c = k c 2 springs s u ↑ -u ↓ 2 cos(θ u -θ s ) = N s k c 2 u ↑ -u ↓ 2 cos 2 (θ u -θ s ) , ( S9
)
where N s is the total number of springs in the system and where the average ⟨•⟩ is performed over all six possible orientations of the coupling springs, namely θ s = iπ/3 with i = 1..6. This averaging of the square cosine yields a result that is independent of θ u . Finally, dividing F c by the total area of the system and noting that there are six coupling springs per hexagon, we find a coupling energy per unit area

f c = 2 √ 3k c 2 u ↑ -u ↓ 2 . ( S10
)
Identifying this expression to the last term of Eq. ( 4) of the main text, we thus find a continuum coupling constant

κ c = 2 √ 3k c . (S11)
In the main text, we use Eqs. (S6) and (S11) to compute values of ν, Γ and ℓ associated with discrete particles and compare the results of our numerical results to continuum predictions in Figs. 3 and4.

S2. CONTINUUM LIMIT COMPUTATIONS

Here we derive the displacement fields of Eqs. ( 9) and ( 12) of the main text from the energy functional presented in Eq. ( 4) of the main text. In Sec. S2 A we derive the general form of the force balance equations associated with the yellow and red sheets. In Sec. S2 B, we solve the force balance equations in a fiber geometry. Then we solve them in a disk geometry in Sec. S2 C. The expressions of the aggregate energies are computed by inserting these results into Eq. ( 4) of the main text and performing the integration. The bulk elastic energy is trivially derived from either one of the resulting expressions by taking the infinite-size limit.

A. Force balance equations

We differentiate the sheet energy density of (S1) with respect to the linearized strain tensor to obtain the constitutive equations of the yellow and red sheets

σ ↑ αβ = λ(∂ γ u ↑ γ -2ϵ)δ αβ + µ(∂ α u ↑ β + ∂ β u ↑ α -2ϵδ αβ ) (S12a) σ ↓ αβ = λ(∂ γ u ↓ γ + 2ϵ)δ αβ + µ(∂ α u ↓ β + ∂ β u ↓ α + 2ϵδ αβ ), (S12b)
where λ = µν/(1ν) is the first Lamé coefficient and where we have used the strain notation of (S2). Differentiating Eq. ( 4) of the main text with respect to the yellow and red sheet displacements u ↑ α and u ↓ α respectively yields the force balance equations for the yellow and red sheets:

∂ β σ ↑ αβ = -κ c (u ↓ α -u ↑ α ) (S13a) ∂ β σ ↓ αβ = -κ c (u ↑ α -u ↓ α ), (S13b)
whose right-hand sides represent the areal densities of external forces exerted by one sheet onto the other through the coupling springs. We parametrize all displacements and stresses by the bulk position vector r. This vector is defined as the position of a point in the bulk state, i.e., in the state characterized by u ↑ = u ↓ = 0. In other words, the actual position of any point of the yellow sheet is given by r + u ↑ (r), and that of a point of the red sheet is given by r + u ↓ (r). In the following, we endeavor to solve the system of equations Eqs. (S12-S13) for the displacement fields u ↑ (r), u ↓ (r) on a two-dimensional domain Ω with the stress-free boundary condition

∀r ∈ ∂Ω n α (r)σ ↑ αβ (r) = n α (r)σ ↓ αβ (r) = 0, (S14)
where n(r) denotes the normal to the domain at a point r of the domain boundary ∂Ω. We use linear elasticity throughout.

B. Fiber

To study the elastic energy of an infinitely long, straight fiber, we write the position vector r = (x, y) in cartesian coordinates. We then solve the force balance equations over the domain (x, y) ∈ Ω = [-W/2, W/2] × R, where W denotes the width of the fiber.

By translational symmetry, the displacement fields gradients ∂ α u ↑/↓ β may only depend on the horizontal coordinate x. this implies that the vertical displacements u ↑ y and u ↓ y are affine functions of the vertical coordinate y. To prevent a divergence in the coupling spring energy, the slopes of these functions must moreover be identical, and we thus write

u ↑ y = u ↓ y = ϕy, ( S15
)
where ϕ is an undetermined constant that cannot depend on x lest diverging strains appear in regions of large y in either or both sheets. The horizontal displacements u ↑ x , u ↓ x must be independent of y for the same reason. Additive constants on the right-hand-side of (S15) can be ignored without loss of generality through suitable choices of the origins of y and u ↑/↓ . As a consequence of (S15), the vertical component of the force exerted by the coupling springs vanishes everywhere.

Inserting these results into Eqs. (S12-S14) yields a system of coupled equations for two functions of one variables, namely u ↑ x (x) and u ↓ x (x):

∂ 2 x u ↑ x = κ c λ + 2µ (u ↑ x -u ↓ x ) (S16a) ∂ 2 x u ↓ x = κ c λ + 2µ (u ↓ x -u ↑ x ), (S16b)
with boundary conditions

∂ x u ↑ x (±W/2) = (λ + µ)2ϵ λ + 2µ - λϕ λ + 2µ (S17a) ∂ x u ↓ x (±W/2) = - (λ + µ)2ϵ λ + 2µ - λϕ λ + 2µ , (S17b)
which is a continuum version of the one-dimensional toy model of the main text, except for the unknown constant ϕ.

To determine the value of ϕ, we note that the fiber as a whole is not subjected to any external force. This implies that its total vertical tension must be constant. Due to the no-stress boundary condition at the y = ±∞ ends of the fiber, this constant is moreover equal to zero:

0 = ∞ -∞ σ ↑ yy (x) + σ ↓ yy (x) dx = 2(λ + 2µ)W ϕ + λU (W/2) -λU (-W/2), (S18)
where we have defined U (x) = u ↑ x (x) + u ↓ x (x) and where the last equality was obtained by using (S13) and performing the integration. We finally combine Eqs. (S16-S17) to obtain a simple differential equation for U (x):

∂ 2 x U = 0 with ∂ x U (±W/2) = - 2λϕ λ + 2µ , (S19)
which implies U (x) = -2λϕx/(λ + 2µ). Inserting this result into (S18) yields ϕ = 0. Inserting the condition ϕ = 0 into the system Eqs. (S16-S17), we find a linear system of differential equations without any unknown parameters. This system thus has a single solution, which can easily be verified to be Eq. ( 9) of the main text.

C. Disk

To study the elastic energy of a disk, we write the position vector r = (r, θ) in polar coordinates. We then solve the force balance equations over the domain (r, θ) ∈ Ω = [0, R] × [0, 2π), where R denotes the radius of the disk.

The rotational invariance of the problem imposes u ↑ θ = -u ↓ θ = 0, and implies that the radial displacement depends only on the radial coordinate. We must thus solve for two scalar functions of one variable, namely u ↑ r (r) and u ↓ r (r). Combining Eqs. (S12-S14), we obtain

∂ 2 r u ↑ r + ∂ r u ↑ r r - u ↑ r r 2 = - κ c λ + 2µ (u ↓ r -u ↑ r ) (S20a) ∂ 2 r u ↓ r + ∂ r u ↓ r r - u ↓ r r 2 = - κ c λ + 2µ (u ↑ r -u ↓ r ) (S20b)
with boundary conditions

∂ r u ↑ r (0) = ∂ r u ↓ r (0) = 0 (S21a) σ ↑ rr (R) = (λ + 2µ)∂ r u ↑ r (R) + λ u ↑ r (r) -ϵ R = 0 (S21b) σ ↓ rr (R) = (λ + 2µ)∂ r u ↓ r (R) + λ u ↓ r (r) + ϵ R = 0. ( S21c
)
By summing these equations two by two, we find a system of equations for U (r) = u ↑ r (r) + u ↓ r (r), namely

∂ 2 r U + U r - U r 2 = 0 with ∂ r U (0) = 0 and (λ + 2µ)∂ r U (R) + λ U (r) R = 0. (S22)
This implies U (r) = 0 and therefore u ↑ r (r) = -u ↓ r (r). Plugging this condition into (S20a) yields

∂ 2 r u ↑ r + ∂ r u ↑ r r - u ↑ r r 2 = 2κ c λ + 2µ u ↑ r , (S23)
which alongside the boundary conditions on u ↑ r comprised in (S21) form a fully specified second-order linear differential equation, for which Eq. ( 12) of the main text is the unique solution.

S4. PROCEDURE TO DRAW RANDOM MATRIX

Here we describe the procedure we use to generate the random instances of the matrix M introduced in Eq. ( 14) of the main text. The form of the energy chose in Eq. ( 14) is very generic, as it boils down to a small-displacement Taylor expansion of any energy function of the 12 vertex coordinates under the constraints of translational and rotational invariance. Here we discuss the way in which we enforce the additional symmetries of the matrix M.

As discussed at the end of the main text, we demand that the elasticity of our particles be three-fold symmetric, i.e., invariant under the permutation:

d 1 → d 2 d 2 → d 3 d 3 → d 1 d 4 → d 5 d 5 → d 6 d 6 → d 4 d 7 → d 8 d 8 → d 9 d 9 → d 7 ,
where the distances d i are defined in Fig. 5(a) of the main text and the permutation above applies simultaneously to the vectors d and d 0 of actual and resting positions. To enforce this condition, we first draw all entries of a 9 × 9 matrix M 0 as independent identically distributed variables from the normal distribution N (0, 1). We then define the 9 × 9 block matrix Ω that enforces the aforementioned permutation as:

Ω =   ω 0 0 0 ω 0 0 0 ω   ,
where the block ω is given by ω

=   0 0 1 1 0 0 0 1 0   . (S24)
We then apply permutations to the symmetry-less matrix M 0 to obtain which has the required three-fold symmetry property.

M 1 = 1 3 M 0 + ΩM 0 Ω -1 + Ω 2 M 0 Ω -2 , (S25)
Our second and last requirement for our elasticity matrix is that it be semi-positive, which prevents the ground state of our individual particles from being mechanically unstable. We enforce this condition by defining

M = M 1 M T 1 , (S26)
where T denotes the usual matrix transposition. By combining Eqs.

[S25] and [S26] and realizing that Ω 3 = I it is easy to show that the energy of Eq. ( 14) of the main text then satisfies the three-fold symmetry condition e d (Ωd) = e d (d).

All elasticity matrices used in Figs. 5 and6 of the main text are obtained through the procedure described here.

S5. SELF-LIMITED AGGREGATE SIZES SCALE LIKE THE BOUNDARY LAYER THICKNESS ℓ IN RANDOM PARTICLES

Randomly drawn matrices M display heterogeneous elastic properties, which result in a wide distribution of aggregate sizes. In Fig. 2 we show that the value of the surface tension Γ/Γ max is not sufficient to accurately predict the size of the aggregates resulting from a matrix M. By contrast, we show in Fig. 5(d) of the main text that rescaling the aggregate sizes by the boundary layer thickness ℓ leads to a collapse of the aggregate sizes, consistent with our continuum theory.

S6. SELECTION OF THE RANDOM PARTICLES USED IN THE PHASE DIAGRAMS

To generate the random particle aggregation diagrams of Fig. 6 of the main text, we first draw 3 × 10 7 random matrices the boundary layer length of the associated particles. This large sample size is required to obtain a sufficient number of particles with relatively high values of ℓ, as is apparent from the fast decay of the probability density of Fig. 3(a) as ℓ increases.

We next divide the range of accessible boundary layer lengths into three intervals: ℓ ∈ [0, 1), ℓ ∈ [1, 1.5) and ℓ ∈ [1.5, ∞). We randomly select batches of 10 6 particles from the first two intervals and use the whole third batch, which contains only 1.25 × 10 5 particles. We then compute the Poisson ratio (see Sec. S7 for the procedure) for all particles in the three batches and further select a few thousand particles from each to obtain quasi-uniform distributions of Poisson ratios as represented in Fig. 3(b). Finally, for each particle, we construct an aggregation diagram by numerically determining the best aggregate upon varying the surface tension. The outcome of this procedure is Fig. 6 of the main text.

S7. PREDICTING THE BOUNDARY LAYER THICKNESS FROM THE ELASTIC PROPERTIES OF RANDOM PARTICLES

In our continuum model, the thickness of the boundary layer that marks the transition between strongly constrained bulk particles and relatively unconstrained aggregate-edge particles is directly tied to the ease with which the two subtriangles that constitute the particles can be shifted relative to each other. Here we tentatively apply a similar reasoning to random particles, and construct the ultimately successful estimate ℓ pred of the boundary layer thickness presented in Fig. 5(c) of the main text. In our continuum model the boundary layer thickness ℓ is constructed from a ratio of elastic moduli, namely

ℓ 2 = λ + 2µ 2κ c = K (1 + ν)κ c . (S27)
where κ c denotes the inter-sheet coupling constant, K = λ + µ is the intra-sheet bulk modulus and ν = λ/(λ + 2µ) is the intra-sheet Poisson ratio. These parameters are not rigorously well-defined in a bulk aggregate of random particles characterized by an elasticity matrix M [Eq. ( 14) of the main text], which does not in general exactly map onto the continuum energy of Eq. ( 4) of the main text. To nonetheless derive our estimate ℓ pred , here we set out to compute proxies for each of these three parameters. We base our procedure on the computation of pseudo-moduli associated with specific deformations illustrated in Fig. 4. We thus define the pseudo-modulus K associated with a deformation vector δd = dd 0 as

K[δd] = δd T • M • δd A hex 0 , (S28)
where A hex 0 = √ 3/2 is the resting area of a hexagonal particle to zeroth order in ϵ. In the following, we choose the normalization of δd so that the definition of (S28) coincides with the moduli discussed in Eq. ( 4) of the main text when applied to our simple particle model (Fig. 2 of the main text). In that specific case, the deformations δd used below are eigenvectors of M.

We first estimate the sheet-coupling modulus κ c as the pseudo-modulus associated with a relative displacement between the two subtriangles of a particle. We picture two possible protocols for such a displacement in Fig. 4(a) and (b), namely a horizontal or a vertical shift of the two triangles, although any intermediate direction is also allowed. Due to the three-fold symmetry of the particle, all these shifting protocols are associated with the same modulus. Using the deformation mode of Fig. 4(a), we write

κ pred c = K 0, 0, 0, 0, 0, 0, - √ 3/2, 0, √ 3/2 . (S29)
We next estimate the intra-sheet bulk modulus K, which characterises the stiffness of each sheet with respect to an isotropic expansion or compression, as the pseudo-modulus associated with the expansion-compression deformation illustrated in Fig. 4(c), namely

K pred = K -2 -3/2 , -2 -3/2 , -2 -3/2 , 2 -3/2 , 2 -3/2 , 2 -3/2 , 0, 0, 0 . (S30)
Our third step is to estimate the Poisson ratio as the average (b) Another displacement protocol equivalent to that of the previous panel. (c) Simultaneous isotropic bulk deformation of the two triangles used to estimate the sheet bulk modulus K pred . In the case of the simple particles of Fig. 2 of the main text, this mixture of expansion and compression imposes bulk deformations on both sub-triangles while leaving the coupling springs lengths unchanged to lowest order in deformation. We thus use it to extract the value of the bulk modulus independently from the coupling modulus in (S30). (d) Bulk and shear deformation of the yellow triangle used to compute ν ↑ through K ↑ and µ ↑ in (S32). The grey nodes are assumed to be free to move in such a fashion that the length of the grey segments remains unchanged during the deformation. This is equivalent to decoupling the black triangle from the rest of the particle, as mentioned in the main text. (e) Bulk and shear deformation of the red triangle used to compute ν ↓ through K ↓ and µ ↓ in (S32). The meaning of the grey nodes is the same as in the previous panel.

ν pred = ν ↑ + ν ↓ 2 (
of the individual pseudo-Poisson ratios ν ↑ and ν ↓ of the yellow and red sublattices (and thus of the yellow and red sheets in the continuum limit). We estimate each of these sheet-specific Poisson ratios as

ν ↑/↓ = K ↑/↓ -µ ↑/↓ K ↑/↓ + µ ↑/↓ , (S32)
where we define the sheet-specific bulk and shear pseudo-moduli through the four deformations illustrated in Fig. 4(de), namely K ↑ = K [(1/2, 1/2, 1/2, 0, 0, 0, 0, 0, 0)] (S33a)

µ ↑ = K 0, - √ 3/4, √ 3
/4, 0, 0, 0, 0, 0, 0

K ↓ = K [(0, 0, 0, 1/2, 1/2, 1/2, 0, 0, 0)] (S33c)

µ ↓ = K 0, 0, 0, √ 3/4, 0, -√ 3/4, 0, 0, 0 .

We finally combine Eqs. (S29-S31) by inserting them into (S27), which yields the values of ℓ pred displayed in Fig. 4(c) of the main text. We also use the pseudo-Poisson ratio defined in (S31) as the vertical coordinate of the phase diagrams of Fig. 6 of the main text.

S8. UNSUCCESSFUL ALTERNATIVE PREDICTOR OF THE BOUNDARY LAYER THICKNESS IN RANDOM PARTICLES

The successful method described in Sec. S7 is only one of many possible extensions of our continuum theory to random particles. Here we discuss a different, unsuccessful approach and draw conclusions from its failure. The foundation of our continuum approach is the existence of an emergent length scale in our deterministic 2D particle model that diverges in the limit where one elastic constant of the particles (k c ) becomes much smaller than another (k). Mechanistically, this large mismatch in elastic constants means that the restoring forces from the softer deformation mode of the particle slowly accumulate from the edge of the aggregate to the bulk as discussed in the main text for the one-dimensional model of Fig. 1. Mathematically, any deformation of the particle may be decomposed into a linear combination of the eigenvectors of the matrix M defined in Eq. ( 14) of the main text. These deformation eigenmodes do not couple to each other in an isolated particle, and each has its own stiffness associated with the corresponding eigenvalue of M. In our deterministic 2D particles, the length ℓ is inversely proportional to the square root of the smallest of these eigenvalues.

This leads us to hypothesize that any soft mode of deformation of the particle may be able to play the same role that the shifting mode of Fig. 4(a-b) plays in our deterministic 2D model: to produce slowly accumulating stresses that take the particle from a more relaxed edge configuration to the bulk configuration. If several soft modes exist within the particle, the softer one should correspond to the thickest boundary layer and thus should dominate the decay of the elastic deformation far enough from the aggregate edge. This reasoning thus suggests the following proxy for the boundary layer thickness:

ℓ eigen = 1 min i λ i , (S34)
where the λ i denote the eigenvalues of M. We test the accuracy of this predictor in Fig. 5 using the same plotting convention as in Fig. 5(c) of the main text, and find that it does not significantly correlate with ℓ. As detailed in the discussion section of the main text, we conclude that not all soft modes of our particles are compatible with the mechanism of stress accumulation over large length scales described above.

S9. TRIANGULAR PARTICLES

To demonstrate the robustness of the results of the main text to a change in particle design, here we introduce a model of triangular deterministic particles. Similar to the main text, the model particles comprise two triangles made of hard k springs and a set of softer (six in this case) k c coupling springs. The particle design is illustrated in Fig. 6(a), and aggregates thereof are shown in Fig. 6(b) and (c). This design gives rise to the same continuum theory as the model of the main text, and we conduct comparisons of pre-made discrete aggregate structures as well as Monte-Carlo simulations using the same methodology as in the main text. These results are shown in Fig. 6(d). The aggregate designs used for the former type of analysis are shown in Fig. 6(e). By contrast with the model described in the main text, the bulk with hole is never advantageous because sticking two triangular particles together has a non-zero elastic energy cost. As a result, a gas of isolated single particles is favored at very low tensions. In addition, unlike in Fig. 3(b) of the main text, the frontier between fiber and bulk is perfectly vertical. This further confirms that the extended region of stability of the fibers observed for the hexagonal model of the main text is due to the curvature of the fibers. Apart from these nuances, the results obtained with this model are very similar to those of the deterministic 2D hexagonal particle model detailed in the main text, thus confirming the robustness of our continuum approach.

[1] H. S. Seung and D. R. Nelson, Defects in flexible membranes with crystalline order, Physical Review A 38, 1005 (1988).

FIG. 2 .

 2 FIG. 2. 2D model of ill-fitting self-assembly. (a) Individual particles have two kind of vertices, which define three interconnected elastic networks. (b) As the particles are aggregated by matching the vertices with the same color, individual particles undergo deformations.

FIG. 3 .

 3 FIG. 3. Large binding energies favor bulk aggregation and incompressible particles tend to form fibrous aggregates. (a) Analytical phase diagram derived from Eqs. (8), (11) and (13). Fiber widths always diverge when approaching the transition to the bulk. Conversely, in disks the radius discontinuously jumps from a finite value to +∞ at the transition. The point where the three phases meet is Γ = 1, ν = 1/2. (b) Phase diagram based on the numerical comparison of the energies of the aggregates shown in Fig. 4(a) for ℓ = 5. The fiber region of the phase diagram is larger than in the continuum model, and the disk radius R again jumps discontinuously at the transition with the bulk. Smaller values of ℓ lead to a very similar phase diagram, albeit with an extended "bulk with holes" region (Supporting Information). Monte-Carlo simulations for the conditions indicated by the small squares are shown as small panels, and are consistent with the phase diagram. The bottom line of snapshots uses 300 particles and a box of 30 × 30 sites, while the others use 200 particles and a box of 60 × 60 sites.

FIG. 5 .

 5 FIG.5. The aggregation behavior of particles with random elasticity is accurately captured by our continuum theory. (a) The shape of a particle is fully characterized by a tuple of nine distances (d1, . . . , d9). (b) To numerically compute the boundary layer thickness ℓ, we simulate a half-plane full of particles (colored here by absolute deformation energy). We then perform an exponential fit of the excess elastic energy relative to the bulk as a function of distance from the edge, shown here as a line on a lin-log plot. Finally we define ℓ as the associated decay length. (c) The random particle boundary layer thickness ℓ measured using this protocol is well predicted by the value ℓ pred inferred from an extrapolation of (10) (red line). (d) Radius and width of the best disk and fibers obtained for random particles. The continuum predictions are shown as red lines (the position of the line does not depend on ν for fibers).

Figure 3

 3 Figure 3 of the main text presents the phase diagram of our non-random discrete particles for boundary layer thicknesses ℓ = +∞ [Fig. 3(a), continuum theory] and ℓ = 5 [Fig. 3(b), numerical procedure described in the main text]. To further illustrate the influence of the boundary layer thickness ℓ, in Fig.1we follow the same procedure used to compute the diagram of Fig.3(b) for a smaller value ℓ = 2.5. Consistent with the results obtained for random particles in Fig.6of the main text, we observe that smaller values of ℓ induce a loss of disks to the benefit of bulks with holes, but no dramatic changes in the fiber and bulk regions of the diagram.

FIG. 1 .

 1 FIG.1. Phase diagram based on the numerical comparison of the energies of the aggregates shown in Fig.4(a) of the main text for ℓ = 2.5. The color code is as in Fig.3of the main text; black: bulk with holes, red: disk, blue: fiber, purple: bulk.

FIG. 2 .

 2 FIG. 2. Equilibrium width W * and radius R * of random particles aggregates in units of number of particles. The data shown is identical to that of Fig. 5(d) of the main text, only without rescaling by the boundary layer thickness ℓ.

FIG. 3 .

 3 FIG. 3. Distribution of random particle properties resulting from the procedure of Sec. S4. (a) Distribution of boundary layer lengths within the initial draw of 3 × 10 7 elasticity matrices M. (b) Distribution of Poisson ratio within each of the three batches of particles. We respectively pick 3400, 3500 and 2400 particles out of the three batches to obtain the quasi-uniform distributions of Poisson ratios materialized by the black lines.

FIG. 5 .

 5 FIG.5. The predictor ℓ eigen of (S34) does not correlate with the measured boundary layer thickness ℓ. Each of the two side panels is a close-up of a region of the central panel, as indicated by the black rectangles.

FIG. 6 .

 6 FIG. 6. Alternative triangle-based model for frustrated self-assembly. (a) Particle design as in Fig. 2(a) of the main text. (b) Example of mechanically equilibrated aggregates as in Fig. 2(b) of the main text. To improve visibility, coupling (blue) springs are not represented. (c) In an infinite (bulk) aggregate, the length of the edge of the two triangles must match. (d) Phase diagram and Monte-Carlo simulation results for ℓ = 2.5 as on the right-hand-side of Fig. 3 of the main text. The bottom line of snapshots contain 300 particles in a box of 30 × 30 lattice sites. The other snapshots contain 200 particles in boxes of size 40 × 40. (e) Examples of aggregate topologies whose energies we compare to compute the phase diagram.

ACKNOWLEDGMENTS

ML was supported by Marie Curie Integration Grant PCIG12-GA-2012-334053, "Investissements d'Avenir" LabEx PALM (ANR-10-LABX-0039-PALM), ANR grants ANR-15-CE13-0004-03, ANR-21-CE11-0004-02, ANR-22-ERCC-0004-01 and ANR-22-CE30-0024-01, as well as ERC Starting Grant 677532. ML's group belongs to the CNRS consortium AQV.

M. Mert Terzi and Hugo Le Roy contributed equally to this work. MMT performed the analytical calculations and HLR conducted the numerical simulations. All authors participated in designing the project and writing the paper.

Supplementary material Collective deformation modes promote fibrous self-assembly in protein-like particles

Hugo Le Roy, 1, 2, * M. Mert Terzi, 1, 3 and Martin Lenz 1, 3, † 1 Université Paris-Saclay, CNRS, LPTMS, 91405, Orsay, France 2 Institute of Physics, École Polytechnique Fédérale de Lausanne-EPFL, 1015 Lausanne, Switzerland 3 PMMH, CNRS, ESPCI Paris, PSL University, Sorbonne Université, Université de Paris, F-75005, Paris, France

S1. MAPPING BETWEEN DISCRETE AND CONTINUOUS PARAMETERS

In this section, we give the relationship between the parameters of our simple (non-random) discrete particle model and the continuum model. The former has three parameters: two spring constants k and k c , and an area stiffness k area . The continuum limit of this model comprises two coupled elastic sheets corresponding to the colors yellow and red in Fig. 2 of the main text, which we respectively denote by the ↑ and ↓ symbols. We represent the elasticity of each sheet by a shear modulus µ and a Poisson ratio ν. The elastic coupling between the sheets is parametrized by the coupling constant κ c . Here we determine µ, ν and κ c in terms of k, k c and k area .

We first map the energy of a single triangular spring network in the discrete particle model onto the energy of a single sheet in the continuum model. The corresponding continuum sheet energy density reads

where u αβ denotes the linearized strain tensor and summation over repeated indices is implied. This strain is expressed with respect to the resting configuration of the sheet of interest, and not with respect to the bulk configuration as in Eq. ( 4) of the main text. The connection between the two conventions is given by the substitutions

ϵδ αβ for the yellow sheet (S2a)

To lowest order in ϵ, the expression of the energy density of (S1) does not depend on whether it is defined as the energy per unit surface of the sheet in its own resting state or in the bulk state defined in the main text. The energy of the discrete triangular spring network is the sum of two parts: the springs and areal energies. The springs part yields the following contribution to the sheet energy density [START_REF] Amos | Arrangement of subunits in flagellar microtubules[END_REF]

Now focusing on the areal stiffness part, we find that its contribution to the total sheet energy reads

where A 0 is the area of a triangle and where we have used the small-deformation approximation of the relative area change as the trace of strain tensor: (A -A 0 )/A 0 ≈ u αα . The corresponding energy density in the continuum limit then reads f area = k area 4 (u αα ) 2 , (S5) * h.leroy@epfl.ch † martin.lenz@universite-paris-saclay.fr arXiv:2308.04698v1 [cond-mat.soft] 9 Aug 2023