Nonequilibrium self-assembly can be found in various biological processes where chemical potential gradients are exploited to steer the system to a desired organized structure with a particular function. Microtubules, for example, are composed of two globular protein subunits, *α*-tubulin and *β*-tubulin, which bind together to form polar dimers that self-assemble a hollow cylinder structure in a process driven by GTPase activity. Inspired by this process, we define a generic self-assembly lattice model containing particles of two subunits, which is driven out-of-equilibrium by a dimer-favoring local driving force. Using Monte Carlo simulations, we characterize the ability of this system to restore pre-encoded target structures as a function of the initial seed size, interaction energy, chemical potential, number of target structures, and strength of the nonequilibrium drive. We demonstrate some intriguing consequences of the drive, such as a smaller critical seed and an improved target assembly stability, compared to the equilibrium scenario. Our results can expand the theoretical basis of nonequilibrium self-assembly and provide deeper understanding of how nonequilibrium driving can overcome equilibrium constraints.

## I. INTRODUCTION

Self-assembly processes play a major role in every living system. Biology has given us numerous examples of many-body self-assembly, including the assembly of the ribosome,^{1} protein folding,^{2} virus capsids,^{3} the nuclear pore complex,^{4} and actin filaments.^{5} Another key example, which is involved in various cellular functions, is microtubules—the most rigid form of cytoskeletal elements, built from 13 parallel filaments forming a tube. These filaments, called *protofilaments*, are chains of dimers, where each dimer is composed of *α*-tubulin and *β*-tubulin protein subunits.^{6}

The astonishing examples of self-assembly nature has accomplished, encourage us to gain deeper understanding of such emergent organization in complex many-body systems. Analyzing generic model systems performing self-assembly can provide insights into biological processes, such as the assembly of protein complexes,^{7} or possible universal aspects of genotype–phenotype maps.^{8} Practically, such minimal models may provide ways for creating desirable molecular structures in practice, e.g., coating particles with certain patches for assembling rings, pyramids, or tetrahedra.^{9} In addition, simulation of a relatively simplified model of the self-assembly of DNA origami can be performed to study various thermodynamical aspects of different origami designs.^{10}

Minimal models have also been used in the investigation of active matter: nonequilibrium locally driven multicomponent systems, in which interacting particles perform collective directed motion.^{11–14} Simulating a minimal model describing such collectively moving particles has shown a phase transition between a disordered and collective movement, estimated its critical exponents, and led to an insight into the applicability of the theoretical framework of critical phenomena to such a nonequilibrium process.^{15} A different schematic model was used to obtain hydrodynamic equations for density and velocity fields of the particles, which are consistent with previous phenomenological results.^{16}

Consequently, extensive research has been performed regarding self-assembly systems and their ability to restore multiple target structures as a function of temperature,^{17} interaction topology,^{18} interaction strength,^{7,19–21} cooperativity,^{22} chemical potential,^{7,19,20} or the concentrations of components.^{19,20,23}

In many cases, biological processes occur far-from-equilibrium, rendering the principle of free-energy minimization irrelevant. Such processes, for example, govern the dynamic instability of microtubules.^{24} The tubulin dimers forming the microtubules bind guanosine-5′-triphosphate (GTP), where the *β*-tubulins hydrolyze GTP to guanosine 5′-diphosphate (GDP), destabilizing the assembly. Since the binding of more dimers to the protofilament catalyzes the GTP hydrolysis, the structure becomes less stable up to a *catastrophe*, where the microtubule falls apart and begins to shrink. In turn, the disassembled free dimers return to the GTP-state, promoting assembly and the *rescue* of the microtubule.^{6} This is a nonequilibrium process fueled by the chemical potential difference between GTP and GDP in the cell cytosol and the GTPase activity of the *β*-tubulins in the protofilament.

Many studies have demonstrated experimentally nonequilibrium self-assembly driven by external forces. For example, the driving may be done by chemical fuel,^{25} light,^{26} or fluid flow.^{27,28} Moreover, experiments have demonstrated the assembly of crystals in rotating magnetic fields,^{29} GTP-driven assembly of microtubules,^{30} and pH or temperature triggered self-assembly of dynamic covalent surfactants.^{31}

From the theoretical perspective, several studies have focused on nonequilibrium self-assembly.^{32} A toy lattice model with a local nonequilibrium driving force has shown improvements in both stability and assembly speed compared to the equilibrium dynamics.^{33} Other systems, driven out-of-equilibrium by a global chemical potential driving force, were found to obey a set of thermodynamic bounds.^{34} In addition, nonequilibrium driving by time-oscillatory inter-particle potentials was shown to result in nontrivial dissipative structures,^{35} and a simple protocol of time-dependent interactions was shown to improve self-assembly in a model colloidal system.^{36}

Inspired by the structure of microtubules, we present a generic toy model of a system having two-subunits coupled to heat and particle baths, which assembles pre-encoded target structures, driven out-of-equilibrium by a local driving force, which favors dimer formation. We analyze the properties of such a system for different interaction energies, chemical potential values, number of encoded targets, seed sizes, and nonequilibrium driving strengths and show several improvements in the ability of the system to restore pre-encoded target structures, including faster assembly, and a smaller critical seed. We also discuss the limitation of the proposed nonequilibrium drive, which is a higher rate of undesired targets assembly.

## II. MODEL PROPERTIES

Our model consists of *N* species of particles: {*i* = 1, 2, …, *N*}, where we assume that *N* is even and a square number. The species are divided into two groups: The first *N*/2 species are called *α**-subunits* and the rest are called *β**-subunits*. A particle of species *i* is denoted as *s*_{i}. The particles are placed on a 2D square lattice consisting of *N*_{0} sites (*N*_{0} is a square number), where each site can be either occupied by a single particle or empty. Each particle can interact with up to *z* other particles (nearest neighbors, denoted as n.n.), where the directionality of the bond of *s*_{i} with *s*_{j} is defined as *δ* = 1, 2, …, *z*. The system operates in the grand-canonical ensemble, with temperature *T*, and all species share the same chemical potential *μ*. In each simulation, there are *m* target structures, each one consists of *M* particles, where we choose *N*_{0} = 4*M*. Throughout the simulations in this work, we used *N*_{0} = 400, which is also the maximal number of particles that can be simultaneously present along the realization.

For our simulations, we take *M* = *N* and *z* = 4, which means that *δ* = 1, 2, 3, 4 stands for the directions *up, right, down*, and *left*, respectively. For simplicity, we take *k*_{B}*T* ≡ 1. This means that we measure *E* and *μ* in units of *k*_{B}*T*; thus, changing the temperature is equivalent to multiplying *E* and *μ* by a constant, while *μ*/*E* remains constant. We also note that the proposed model is described in the grand-canonical ensemble, namely, the lattice is coupled to both a particle-bath and a heat-bath, with which it can exchange particles and heat. Therefore, the chemical potential of each species is related to the concentration of this species in the particle-bath.

### A. Target structures

We define a *target chain* as a row or a column of $M$ particles of different species, where no two adjacent particles are of the same subunit. For example, for *M* = *N* = 16, a possible target chain is shown in Fig. 1(a).

Each target structure is a $M\xd7M$ square, made of $M$ parallel target chains, where each chain contains different particles. An example of a possible target structure, for the same set of {*M*, *N*} values, is shown in Fig. 1(b).

### B. Energy

In this model, all particle species have zero internal energy. The interaction energy between two adjacent particles *s*_{i}, *s*_{j} where the directionality of *s*_{j} with respect to *s*_{i} is *δ* is given by

Examples for interaction energies between pairs of particles are depicted in Fig. 1(c).

### C. Simulations

We run Monte Carlo (MC) simulations with the Metropolis algorithm.^{37} In each simulation, the initial state of the lattice is the *S* × *S* central region of one of the target structures. This central part is a *seed* of length *S*, and the respective target structure is denoted as the *seeded* or *desired* target. An example for a seed is shown in Fig. 1(d). In each step, one of the *N*_{0} sites of the lattice is chosen randomly. Then, a species from the set {0, 1, 2, …, *N*} is chosen randomly, where “0” stands for “unoccupied.”

At equilibrium, the site accepts the new species with probability *p*_{eq},

where Δ*U* is the difference in energy *U* between the proposed state and the current state, and *U* is defined as

where *N*_{oc} is the current number of occupied sites and the summation is over all pairs of nearest neighbor particles *s*_{i}, *s*_{j} on the lattice.

Now that we have established the model properties at equilibrium, a nonequilibrium local driving force is added.

### D. Driving force

The local driving force is included through a modulation of the acceptance probability of a Monte Carlo step,^{33} *p*_{eq}. The essence of the proposed drive is increasing the acceptance probability of steps which promote the assembly of chains and thus breaking detailed balance. Formally, the driving force encourages the formation of the following structure: two neighbors of a particle along the same axis are of the same subunit, which is the opposite to the particle between them, namely, *α*–*β*–*α* or *β*–*α*–*β* triplets. We define such a structure as an *intermittent subunit triplet*. The new acceptance probability is

where *D* = Δ*λ* if the proposed step completes an intermittent subunit triplet, *D* = −Δ*λ* if the proposed step destroys an existing intermittent subunit triplet, and otherwise *D* = 0. Δ*λ* is, therefore, the drive strength.

Whenever multiple triplets are affected by a proposed step, the value of *D* is determined by a majority voting of *D* values from all triplets involved, where tie means *D* = 0. Formally, we take $D=\Delta \lambda \xd7sign\u2211iDi$, summing over all triplets affected by the proposed move. We note that such a definition accounts for a case where the mechanism creating the local driving force reaches a saturation value when the number of intermittent subunit triplets following a certain move is larger than the current number of intermittent subunit triplets by 1. This way, the saturation effect is realized with reasonable probability, and a driving force with a saturating nature is being modeled. Further elaboration can be found in the supplementary material, Sec. S1. Examples for several scenarios and their respective *D* values are shown in Fig. 2, along with an example for a cycle breaking detailed balance.

## III. QUANTIFYING ASSEMBLY

In order to quantitatively characterize the assembly of a target structure, given a set of pre-defined parameters such as *μ*, *E*, *S*, and *m*, we use the following measures:

*Final distance*(*d*_{f}): distance from the desired target structure, averaged over the last 5% of the simulation steps (5000 lattice sweeps for typical simulations, where one lattice sweep equals*N*_{0}Monte Carlo steps), where the*distance*(*d*) is defined as the difference between*M*and the number of particles that are in their designated location according to the target structure.*First assembly time*(*T*_{fa}): time, measured in Monte Carlo steps, until the system first assembles a structure close to the desired target structure, within a distance determined by the average root mean square (rms) fluctuations around the target. The motivation for the definition of this*fluctuation deviation*and its derivation can be found in the supplementary material, Sec. S3.*Time spent in target*(*T*_{t}): the total time the system is within the rms fluctuation distance around the target.*Chimeric value*(*CV*): the number of particles in the biggest structure, which do not belong to the desired target. The chimeric value is averaged over all the simulation steps. This parameter quantifies the size of chimeric structures bound to the growing seed, where lower CV means a more selective assembly of the desired target.*Characteristic growth time*(*T*_{g}): time until the system first assembles a structure of size*M*. Note that it does not necessarily have to be the target structure, but rather any structure of size*M*. This parameter quantifies the inverse of the aggregate growth rate.

An example for a distance-over-time plot is presented in Fig. 3, where some of the above measures are denoted. We note that Fig. 3 shows the distance from the desired target and thus only partly reveals the simulation dynamics after assembling this target (*T* > *T*_{fa}). More elaboration on post-assembly dynamics can be found in the supplementary material, Sec. S4.

## IV. FINE TUNING OF THE SYSTEM

We start by exploring the {*μ*, *E*} parameter space, where *E* is the binding energy and *μ* is the chemical potential of all of the species, to find the phase space regime in which the system can retrieve the stored target structures.

Qualitatively, too low chemical potential results in a high rate of particles entering the lattice, whether they fit the desired target structure or not. Conversely, too high chemical potential results in a high rate of particles leaving the lattice, hindering the formation of large structures. Similarly, if the binding energy is too weak, there is no significant energetic preference to the target structure compared to a random organization of the particles, whereas too high binding energy leads to kinetic traps along the assembly.^{33}

Roughly, the optimal region in phase space can be narrowed down by simple calculations. Let us denote the microstate *s*_{t}, in which the lattice solely contains a target structure of *M* particles, and define $l\u2261M$. For *s*_{t}, there are *l*^{2} particles on the lattice with 2*l*(*l* − 1) strong bonds. Therefore, according to Eq. (3), the total energy *U* of the microstate *s*_{t} is given by *U*_{t} = −2*l*(*l* − 1)*E* + *l*^{2}*μ*. The total energy *U* of a microstate *s*_{e} corresponding to an empty lattice is 0. In order to achieve assembly, we first choose parameters such that the equilibrium probability of the system being at the microstate *s*_{t} is larger than the probability of an empty lattice, *s*_{e}. Since the equilibrium probability of a certain microstate is proportional to its Boltzmann factor *e*^{−U}, we get *U*_{t} < 0, which gives

Substituting the value *l* = 10, which is used later in the simulations, we get *μ*/*E* < 1.8.

For estimating optimal values for the binding energy *E*, we seek for a condition that differentiates between desired and undesired MC moves and ensures that the undesired ones are accepted with small enough probability during each simulation step. The condition we choose states that removal of the seed’s least bounded particle is more likely to be rejected than accepted. Given a square-shape seed, the most loosely bound particles are the ones in its corners. The probability of accepting an undesired step that removes a corner particle is given by *p*_{r} = exp{ −2*E* + *μ*}. In order for this unwanted move to have a better than an even chance to be rejected, we set *p*_{r} < 1/2. This inequality, together with Eq. (5), determine *E* as follows:

Full derivation of Eq. (6) can be found in the supplementary material, Sec. S5. Substituting the value of *l* = 10, condition *E* > 5 ln(2) ≈ 3.47 is obtained.

As for the driving force value Δ*λ*, we use the following observation: At equilibrium, the undesired step with the highest acceptance probability that is still smaller than unity is the removal of a particle while breaking two strong bonds (e.g., removing a corner particle from the seed). Its corresponding acceptance probability is *p*_{r} = exp{ −2*E* + *μ*}. While applying an encouraging driving force (*D* = Δ*λ*) in this situation is impossible, choosing *e*^{Δλ}*p*_{r} < 1 ensures that the driving force would not be able to make any undesired step, a one which is accepted with probability 1. Namely, we set

To test the above statements and find the optimal region in phase space for the assembly, we ran the simulations and calculated the final distance of the system from the desired target structure, with *N* = *M* = 100, *m* = 6 targets, a seed length of *S* = 7, and a lattice size of *N*_{0} = 400. The energy *E* was varied in the range [5, 12] and the ratio *μ*/*E* was varied in the range [1.2, 2], which roughly correspond to the estimated optimal region in parameter space, and the system was examined both at equilibrium, Δ*λ* = 0, and out-of-equilibrium, taking Δ*λ* = 1.

These simulations yielded a set of (*μ*/*E*, *E*, *d*_{f}) points, for both the equilibrium and nonequilibrium cases, which can be presented as a heat-map of the *d*_{f} values as a function of *μ*/*E* and *E* values, shown in Figs. 4(a) and 4(b) following a Gaussian interpolation^{38} (see the supplementary material, Sec. S6). The interpolation was applied in order to improve the resolution of the obtained heat-map.

From these results, we recognize three regimes in (*μ*/*E*, *E*) space [Figs. 4(a) and 4(b)]. In order to inspect these regimes more carefully, the state of the lattice by the end of the simulation was captured, at both equilibrium and nonequilibrium conditions, for two regions in (*μ*/*E*, *E*) space—at the left boundary of regime III [Figs. 4(c) and 4(d)] and at the right boundary of regime III [Figs. 4(e) and 4(f)].

In regime I, the system partly assembles the desired target structure as the distance values are larger than twice the fluctuation deviation, in agreement with the qualitative prediction of a *μ* value too low. In regime II, the initial seed dissolves and all the particles leave the lattice, as reflected in a distance value close to *M* = 100. Such a distance does not necessarily indicate an empty lattice, but here, the final lattice state was examined and indeed found to be empty. This scenario is in agreement with both the qualitative and quantitative estimations since in the equilibrium heat-map, this regime is bounded by *μ*/*E* ≳ 1.8. In regime III, defined as the region where the distance from the target is lower than two fluctuation deviations, a compromise between the limitations of the two other regimes is found. Hence, it is considered as the optimal regime in parameter space for the self-assembly and will be used for later simulations.

These regimes I, II, and III are found in both equilibrium and nonequilibrium conditions, where in the nonequilibrium case, regime II is smaller and regime III remains approximately the same size. This is a positive effect of the driving force as the initial seed does not dissolve for a wider range of {*μ*, *E*} parameters.

This behavior may be explained by a closer inspection of the influence of the driving force on the lattice growth and the lattice shrinkage. The driving force may or may not encourage a move of a particle entering an empty site but would never discourage it as such a move may or may not form a new intermittent subunit triplet but cannot destroy one. In a similar manner, the driving force would never encourage a move where a particle is removed from the lattice since such a move would never increase the total number of intermittent subunit triplets and may discourage it (when a removal of a particle results in a destruction of an intermittent subunit triplet). Consequently, for high *μ*/*E*, where lattice shrinkage is preferred over growth, the drive cannot encourage shrinkage and may possibly discourage it. This way, the driving force restrains the extensive shrinkage, which characterizes this region in (*μ*/*E*, *E*) space, as can be seen in Figs. 4(c) and 4(d). This effect enables the system to assemble the target structure for higher *μ*/*E* values, which is seen on the heat-maps as regime II becomes smaller and regime III stretches to higher *μ*/*E* values. Conversely, for low *μ*/*E*, where lattice growth is preferred over shrinkage, the drive cannot discourage growth but may encourage unwanted growth (e.g., insertion of a particle to a site where it forms an intermittent subunit triplet but does not correspond to the desired target structure). Thus, the driving force further fuels the extensive growth of undesired structures, which characterizes this region in (*μ*/*E*, *E*) space, as shown in Figs. 4(e) and 4(f). This effect allows for successful nonequilibrium assembly only for higher *μ*/*E* values. The two effects mentioned above make nonequilibrium assembly possible for higher *μ*/*E* values compared to the equilibrium case, while preventing it for lower *μ*/*E* values, resulting in regime III being shifted to the right while retaining its size in the driven scenario presented in Fig. 4(b).

## V. NUMBER OF PRE-ENCODED TARGET STRUCTURES

In order to challenge the capacity of the system at equilibrium and nonequilibrium conditions, we ran the simulations with an increasing number of targets and characterized the time spent in target, *T*_{t}, and the first assembly time, *T*_{fa} (Fig. 5). Specifically, we ran simulations with 10^{5} lattice sweeps (4 × 10^{7} MC steps), with the parameter set {*m* = 2, 3, …, 9, *N* = *M* = 100, *N*_{0} = 400, *E* = 8, *μ*/*E* = 1.78, *S* = 5, Δ*λ* ∈ [0, 1.7]} chosen from the optimal range within regime III.

For *m* > 6, the system does not assemble the desired target structure with or without the driving force, as evident from the values of *T*_{t} = 0 [Fig. 5(a)] and *T*_{fa} > *T*_{tot} [Fig. 5(b)], indicating a storage capacity of *m*_{max} = 6 for the parameters chosen.

The trends of *T*_{t} and *T*_{fa} as a function of the number of encoded targets, showing optimal behavior for intermediate *m* values, can be rationalized through a trade-off between the probability of forming chimeric structures and the mobility of the system. For low *m* values, the number of structures which can form with considerable probability (i.e., parts of target structures) is low, resulting in slow dynamics of the system and slower assembly [Fig. 5(b)]. This, in turn, limits the time the system can spend near the target structure [Fig. 5(a)]. Since the desired result is specifically assembling the seeded target structure, the number of desired structures (e.g., desired pairs and triplets) remains constant as *m* is varied. However, for larger *m*, there is an increasing number of possible structures which are energetically favorable, yet do not belong to the desired target, i.e., chimeric structures, resulting in kinetic traps, which lead to lower assembly quality. When the rate of the formation of chimeric structures is too high, the desired structure cannot be assembled during the simulation.

Although the driving force does not increase storage capacity, as *m*_{max} = 6 also for nonequilibrium conditions, it does contribute to the quality of the assembly for *m* ≤ 6 pre-encoded target structures, which become both more stable (higher *T*_{t}) and assembled faster (lower *T*_{fa}), as seen from Figs. 5(a) and 5(b), respectively.

Next, we chose *m* = 3 to study the effect of the drive on the system. Compared to the equilibrium case, Δ*λ* = 0, nonequilibrium driving both increases the target stability and accelerates the assembly, as evident from the higher values of the time spent in the desired target structure *T*_{t} [Fig. 5(c)] and the lower values of the first assembly time *T*_{fa} [Fig. 5(d)], respectively. However, both *T*_{t} and *T*_{fa} do not vary monotonically by increasing the drive Δ*λ*. Specifically, there exists a value of drive that optimizes both the stability and the assembly speed, Δ*λ* = 0.7, for the given simulation parameter set. This non-monotonic behavior is a consequence of a nonequilibrium trade-off, which is discussed in Sec. VII.

## VI. SUB-CRITICAL SEED

Based on nucleation theory, Murugan *et al.*^{19} calculated the size of a “critical seed” *S*_{c} × *S*_{c} for equilibrium self-assembly, below which the system dissolves,

In order to test this constraint for the nonequilibrium system, we ran simulations with various seed sizes and found the final distance of the system from the desired target structure (Fig. 6). The simulations were run for 10^{5} lattice sweeps (4 × 10^{7} MC steps), with the parameter set {*m* = 6, *N* = *M* = 100, *N*_{0} = 400, *E* = 9, *μ*/*E* = 1.78, *S* = 2, 3, …, 9, Δ*λ* ∈ [0, 1.7]}. Similar to the previous set of simulations, the values of *E* and *μ*/*E* were chosen from the optimal range within regime III. Substituting the values of *μ*/*E* = 1.78, in Eq. (8), we get *S*_{c} = 4.54.

At equilibrium, the final distance from the target is close to 100, the target size, for *S* ≤ 4, meaning that the initial seed is dissolved. For *S* ≥ 5, however, assembly is achieved, in agreement with the theoretical prediction of the critical seed size. In contrast, when the system is driven away from equilibrium, the equilibrium critical seed constraint no longer holds, and the system assembles the seeded target structure for sub-critical seed lengths as well. We again notice an optimal value of the drive Δ*λ* = 1.2, which yields the minimal distance from the target. Again, the existence of an optimal drive value stems from a nonequilibrium trade-off, which is treated in Sec. VII.

For the sake of completeness, we briefly discuss the situation in which no seed is introduced (*S* = 0), and the initial state is an empty lattice. In an empty lattice, a move of a particle entering the lattice is accepted with probability of *p*_{init} = *e*^{−μ}. In order for such a move to happen with a reasonable probability of, for instance, 1%, we should set *μ* ≲ 4.7. In such a case, Fig. 4(a) shows that assembly is unlikely to happen. Namely, for *μ* values low enough, a structure would form (since *p*_{init} is large enough and the rate of particles entering the lattice is high), but its distance from a target structure would be high, and for higher *μ* values, no structure would form. We note that the local driving force does not influence *p*_{init}, as no triplet exists nor can be formed in the initial state, so high-*μ* scenario stays the same also for the nonequilibrium case. For low *μ* values, the lattice grows rapidly, so driving the system would only increase the error, as explained in Sec. IV. Therefore, the drive does not change the above equilibrium conclusion.

## VII. A NONEQUILIBRIUM TRADE-OFF

Our results show a non-monotonic behavior of the time spent in the target *T*_{t}, the first assembly time *T*_{fa}, and the final distance *d*_{f} as a function of the drive. As the drive value increased, the assembly was first accelerated, the targets became more stable, and the assembled structure was closer to the target structure, up to a certain drive value, above which there was a decrease in assembly speed and target stability and an increase in the distance from the target.

In order to unveil the reason behind this phenomenon, we ran simulations for 10^{5} lattice sweeps (4 × 10^{7} MC steps), with a similar parameter set {*m* = 6, *N* = *M* = 100, *N*_{0} = 400, *E* = 7, *μ*/*E* = 1.78, Δ*λ* ∈ [0, 1.9]}, and examined the characteristic growth time *T*_{g} and chimeric value *CV* (Fig. 7). The system was challenged with a high number of pre-encoded targets, *m* = 6, increasing the possibility of forming chimeric structures and thus yielding more valuable information about *CV*.

Indeed, the aggregate growth process is accelerated when increasing the drive, up to certain drive value, Δ*λ* = 1.4, above which the growth rate does not improve anymore. However, the chimeric value is higher for increasing drive values, meaning that in addition to the assembly of the desired target, other chimeras are being assembled as well.

We thus recognize a trade-off between the accelerated assembly and the reduced chimeric value. On the one hand, a stronger drive increases the formation rate of chain-like structures necessary for building the target, thereby accelerating the assembly process. On the other hand, too strong of a drive encourages chain-like structures, which are not necessarily of the seeded target, to bind to the growing seed, reflected in the chimeric value. This trade-off governs the non-monotonic behavior of the time spent in the target, the first assembly time, and the final distance from the desired target as a function of the drive, as shown in Figs. 5 and 6.

## VIII. DISCUSSION

In this work, we have described a toy model for a dimer-based self-assembly system. We examined its ability to assemble pre-encoded target structures at both equilibrium and nonequilibrium regimes. The nonequilibrium drive used in this work increased the acceptance probability of chain-forming moves while decreasing the probability of chain-destroying ones.

Initial exploration of the parameter space yielded three different regimes corresponding to a partial assembly of the seeded target (regime I), dissolution of the seed (regime II), and a successful assembly of the desired target (regime III) as a function of the binding energy *E* and the ratio between the chemical potential and the binding energy *μ*/*E*.

Next, we tested the effect of the nonequilibrium local drive on the assembly quality measured by the time spent in the target, which corresponds to target stability, and the first assembly time, which corresponds to the speed of assembly, and demonstrated an improved assembly compared to the equilibrium case. This improvement, however, does not increase monotonically for larger drive values but rather shows an optimal behavior for intermediate values. The underlying mechanism behind this phenomenon is an inherent nonequilibrium trade-off between the assembly speed and the exclusiveness of the assembled structure, where too fast of an assembly is accompanied by chimeras.^{19} This is an example of a concept present in various biological systems, where increased dissipation, or irreversibility, does not necessarily improve life-supporting properties.^{39}

Perhaps the most striking property of the local drive shown in this work is its ability to promote target assembly in the sub-critical seed region, overcoming the equilibrium constraint. This unveils the power of such nonequilibrium driving force, allowing the system to assemble a target structure in a region in parameter space, which is inaccessible at equilibrium.

This decrease in the critical seed size explains the increase in the rate of chimeric structures compared to the equilibrium scenario. Small seeds of undesired targets are constantly forming, and due to the small critical seed, they grow and create larger chimeras.

While we were inspired by general aspects of the microtubule structure, our model is a generic toy model, aimed at describing a wide range of dimer-based physical systems. For this reason, it is essential to understand the range of systems for which our results and conclusions are applicable. In the following paragraphs, we will examine the compatibility of the proposed model to real-life systems.

For instance, the concept of programming inter-particle interactions in order to make some structures local minima in the free-energy landscape is analogous to DNA origami, where different DNA strands are chosen such that a certain folding of the scaffold strand is achieved through relaxation of the strands mixture toward equilibrium.^{40} Furthermore, our proposed interaction energy landscape, having sharp local minima for certain structures and none for others, can be found in the modeling of the process of protein binding to DNA. In this process, a protein slides along the DNA, searching for its target site, where it binds the DNA. The interaction energy between the protein and the DNA is spatially (sequence) dependent, and it was shown that an adequate model for such energy landscape involves the protein experiencing, for most of the time, a Gaussian energy-sequence landscape with relatively small variance.^{41,42} This follows the same logic of the interaction-energy matrix of our proposed model, accounting for a −*E* interaction energy between pairs of particles from a certain target structure and no interaction otherwise. In addition, dependence of certain assembly probabilities on their chain structure, which appears in our model driving force, exists in microtubules. In microtubules, a dimer at the end of the filament is more likely to move into a specific, confrontationally different, state when another dimer binds the end of the microtubule. This confrontational change, involving hydrolysis of GTP to GDP in the dimer, promotes disassembly of the microtubule, as part of the nonequilibrium process we discussed in Sec. I.^{6}

However, there are certain aspects of the proposed model that limit the range of real-life systems for which our findings may be applicable.

First, the model we proposed is a lattice model. Generally speaking, such a model is a powerful tool, which can capture much of the behaviors of a compatible off-lattice model for patchy particles.^{43,44} When comparing a model of a driven fluid, similar to the one suggested here, with its continuous off-lattice counterpart, it has been shown that while both experience an order–disorder phase transition, several aspects, such as the critical exponent or critical temperature, vary non-trivially.^{45} This much richer phase diagram for an off-lattice model shown by Díez-Minguito *et al.* may imply that while the existence of the critical seed and its relative change compared to the non-driven case are likely to be consistent with an off-lattice model, specific values such as the critical seed length should be taken with extra caution.

Second, while our model includes moves of insertion or deletion of a particle to or from a certain lattice site, it does not allow for the movement of a cluster. This might seem incomplete in light of several studies that describe assembly as a process in which smaller clusters of subunits from a target structure (fragments) are formed and interact^{46} and show that assembly of certain structures may follow a mechanism where pairs of multimers hybridize with one another, leading to the creation of larger bonding networks and to extensive bonding of the subunits.^{47} However, the dynamics of the general set of systems described by our model, determined by the external parameters *μ*, *E*, are slightly different. For low *μ* values, the rate of particles entering the lattice is high, hindering the retrieval of the desired target structure. Therefore, the *μ* value must be high enough, resulting in a large critical seed size [according to Eq. (8)]. As stated by Jacobs *et al.*,^{46} the existence of a large critical seed size implies that newly formed small fragments are more likely to dissolve rather than keep growing, and thus, substantial fragments cannot form in empty lattice sites. Therefore, the only stable structure in the system is the growing initial seed. This kind of dynamics, where new structures only bind to the initial seed rather than creating new nucleation seeds, is demonstrated in the images of the system states in Fig. 3. Consequently, if we allowed particles or clusters to move to neighboring sites as part of the simulation dynamics, this would only result in the main cluster (the grown initial seed) moving back and forth along the lattice, not changing the assembly quality parameters.

There are a few other simplifying assumptions in our model: the first is the zero-interaction energy we chose for undesigned interactions. Nevertheless, we have demonstrated that such interaction energy can serve as a model for finite, low enough, weak interaction energies (see the supplementary material, Sec. S7). The second assumption has to do with the fixed finite size of the lattice that was simulated. However, we have shown that different lattice sizes have minimal effect on the behaviors of the model (see the supplementary material, Sec. S8). The third assumption is the lack of rotational moves in the simulation. Still, adding rotations would only make the assembly times longer, leaving the nature of the dynamics unchanged (see the supplementary material, Sec. S9).

Further studies may be done, examining this dimer-based system in an off-lattice model and achieving a better quantitative estimation of the critical values obtained in this work. Future works may also verify the universality of our conclusions to general *n*-mer systems or to larger more complicated dimer-based systems.

Our findings provide insights into dimer-based biomolecular systems operating far-from-equilibrium, which can be used for designing adaptive materials.^{48} The sub-critical seed operation demonstrated in this work may be exploited for compressing information or associative memory models^{20} as the smaller seed can be expanded to the whole target structure by applying the local drive.

## IX. METHODS

## SUPPLEMENTARY MATERIAL

See the supplementary material for detailed information on the simulation dynamics, justifications for assumptions, and data analysis methods, as referenced in the main text.

## ACKNOWLEDGMENTS

Gili Bisker acknowledges the support from the Zuckerman STEM Leadership Program; the Israel Science Foundation (Grant No. 456/18); the Ministry of Science, Technology, and Space, Israel (Grant No. 3-17426); and the Tel Aviv University Center for AI and Data Science (TAD). Adi Ben-Ari and Liron Ben-Ari acknowledge the support from the Excellence Program of the Faculty of Engineering at Tel Aviv University.

This work was supported by the Air Force Office of Scientific Research (AFOSR) under Award No. FA9550-20-1-0426. Research was sponsored by the Army Research Office (ARO) and was accomplished under Grant No. W911NF-21-1-0101. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the Army Research Office or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

A.B. and L.B. equally contributed to this work.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.