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.

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.

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 si. The particles are placed on a 2D square lattice consisting of N0 sites (N0 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 si with sj 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 N0 = 4M. Throughout the simulations in this work, we used N0 = 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 kBT ≡ 1. This means that we measure E and μ in units of kBT; 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.

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).

FIG. 1.

Lattice model for M = N = 16. The blue squares represent α-subunits, the orange squares represent β-subunits, and the numbers represent the different species. (a) A possible target chain. (b) A possible target structure, containing the target chain from (a). The bold line represents the lattice borders. (c) Examples for interaction energies between pairs of particles, assuming that the target structure from (b) is the only one encoded. The interaction energy is non-zero only if the two particles are in the same configuration as in a certain target structure, including the same directionality. (d) A seed of the target structure from (b), which is the initial state of the lattice in our simulations.

FIG. 1.

Lattice model for M = N = 16. The blue squares represent α-subunits, the orange squares represent β-subunits, and the numbers represent the different species. (a) A possible target chain. (b) A possible target structure, containing the target chain from (a). The bold line represents the lattice borders. (c) Examples for interaction energies between pairs of particles, assuming that the target structure from (b) is the only one encoded. The interaction energy is non-zero only if the two particles are in the same configuration as in a certain target structure, including the same directionality. (d) A seed of the target structure from (b), which is the initial state of the lattice in our simulations.

Close modal

Each target structure is a M×M 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).

In this model, all particle species have zero internal energy. The interaction energy between two adjacent particles si, sj where the directionality of sj with respect to si is δ is given by

Uijδ=E,si,sj are n.n. with bond directionality δwith respect to si in any target,0otherwise.
(1)

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

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 N0 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 peq,

peq=min1,eΔUkBT,
(2)

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

U=Uijδ+μNoc,
(3)

where Noc is the current number of occupied sites and the summation is over all pairs of nearest neighbor particles si, sj on the lattice.

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

The local driving force is included through a modulation of the acceptance probability of a Monte Carlo step,33peq. 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

pneq=min1,eΔUDkBT,
(4)

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=Δλ×signiDi, 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.

FIG. 2.

Local driving force. (a)–(c) Examples of situations yielding different D values. The initial state is presented on top, and the proposed step at the bottom. (a) The initial state does not include an intermittent subunit triplet, and the proposed step does; thus, D = Δλ. (b) The initial state includes an intermittent subunit triplet, and the proposed step does not; thus, D = −Δλ. (c) Neither the initial state nor the proposed one includes an intermittent subunit triplet; thus, D = 0. (d) Example of a cycle breaking detailed balance. The acceptance probabilities of the moves appear near the corresponding arrows. The probability of moving from the leftmost state to the rightmost state depends on the path taken: the probability for the upper path is pu = exp{2E − 2μ + Δλ}, while for the lower path pl = exp{2E − 2μ + 2Δλ}, meaning ln(pu/pl) ≠ 0, i.e., broken detailed balance. Here, we assumed that the target structure from Fig. 1(b) is the only one encoded.

FIG. 2.

Local driving force. (a)–(c) Examples of situations yielding different D values. The initial state is presented on top, and the proposed step at the bottom. (a) The initial state does not include an intermittent subunit triplet, and the proposed step does; thus, D = Δλ. (b) The initial state includes an intermittent subunit triplet, and the proposed step does not; thus, D = −Δλ. (c) Neither the initial state nor the proposed one includes an intermittent subunit triplet; thus, D = 0. (d) Example of a cycle breaking detailed balance. The acceptance probabilities of the moves appear near the corresponding arrows. The probability of moving from the leftmost state to the rightmost state depends on the path taken: the probability for the upper path is pu = exp{2E − 2μ + Δλ}, while for the lower path pl = exp{2E − 2μ + 2Δλ}, meaning ln(pu/pl) ≠ 0, i.e., broken detailed balance. Here, we assumed that the target structure from Fig. 1(b) is the only one encoded.

Close modal

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 (df): 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 N0 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 (Tfa): 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 (Tt): 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 (Tg): 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 > Tfa). More elaboration on post-assembly dynamics can be found in the supplementary material, Sec. S4.

FIG. 3.

Distance from the desired target structure as a function of time, measured in Monte Carlo steps. Highlighted on the plot, Tfa is the time until first assembly, Tt is the time spent in the desired target structure, and e is the fluctuation deviation. The images show the state of the lattice at certain timepoints. Structures that belong to the desired target structure are marked in green, and other, undesired, structures are marked in red. The boundaries of the target structure are marked in blue (for more elaboration on the coloring scheme, see the supplementary material, Sec. S2). The chimeric structures, whose size is the chimeric value, are also highlighted. This example dataset was taken from a simulation of 105 lattice sweeps, with parameters {m = 6, M = 100, N0 = 400, E = 9, μ/E = 1.78, S = 5, and Δλ = 1.5}.

FIG. 3.

Distance from the desired target structure as a function of time, measured in Monte Carlo steps. Highlighted on the plot, Tfa is the time until first assembly, Tt is the time spent in the desired target structure, and e is the fluctuation deviation. The images show the state of the lattice at certain timepoints. Structures that belong to the desired target structure are marked in green, and other, undesired, structures are marked in red. The boundaries of the target structure are marked in blue (for more elaboration on the coloring scheme, see the supplementary material, Sec. S2). The chimeric structures, whose size is the chimeric value, are also highlighted. This example dataset was taken from a simulation of 105 lattice sweeps, with parameters {m = 6, M = 100, N0 = 400, E = 9, μ/E = 1.78, S = 5, and Δλ = 1.5}.

Close modal

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 st, in which the lattice solely contains a target structure of M particles, and define lM. For st, there are l2 particles on the lattice with 2l(l − 1) strong bonds. Therefore, according to Eq. (3), the total energy U of the microstate st is given by Ut = −2l(l − 1)E + l2μ. The total energy U of a microstate se 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 st is larger than the probability of an empty lattice, se. Since the equilibrium probability of a certain microstate is proportional to its Boltzmann factor eU, we get Ut < 0, which gives

μE<2l1l.
(5)

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 pr = exp{ −2Eμ}. In order for this unwanted move to have a better than an even chance to be rejected, we set pr < 1/2. This inequality, together with Eq. (5), determine E as follows:

E>l2ln(2).
(6)

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 pr = exp{ −2Eμ}. While applying an encouraging driving force (D = Δλ) in this situation is impossible, choosing eΔλpr < 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

Δλ<2Eμ.
(7)

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 N0 = 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, df) points, for both the equilibrium and nonequilibrium cases, which can be presented as a heat-map of the df values as a function of μ/E and E values, shown in Figs. 4(a) and 4(b) following a Gaussian interpolation38 (see the supplementary material, Sec. S6). The interpolation was applied in order to improve the resolution of the obtained heat-map.

FIG. 4.

(a) and (b) Final distance from the seeded target as a function of μ/E and E at equilibrium (Δλ = 0) (a) and nonequilibrium (Δλ = 1) (b), following Gaussian interpolation. Colors represent final distances in a logarithmic scale. Each data point is averaged over ten realizations (see the supplementary material, Sec. S6). [(c)–(f)] The state of the lattice by the end of the simulation on the boundaries of nonequilibrium regime III, colored as described in Fig. 3. (c) Right boundary (namely, μ/E and E values, which correspond to the right boundary of nonequilibrium regime III) at equilibrium: E = 7, μ/E = 1.9, Δλ = 0. (d) Right boundary out-of-equilibrium: E = 7, μ/E = 1.9, Δλ = 1. (e) Left boundary at equilibrium: E = 8, μ/E = 1.7, Δλ = 0. (f) Left boundary out-of-equilibrium: E = 8, μ/E = 1.7, Δλ = 1.

FIG. 4.

(a) and (b) Final distance from the seeded target as a function of μ/E and E at equilibrium (Δλ = 0) (a) and nonequilibrium (Δλ = 1) (b), following Gaussian interpolation. Colors represent final distances in a logarithmic scale. Each data point is averaged over ten realizations (see the supplementary material, Sec. S6). [(c)–(f)] The state of the lattice by the end of the simulation on the boundaries of nonequilibrium regime III, colored as described in Fig. 3. (c) Right boundary (namely, μ/E and E values, which correspond to the right boundary of nonequilibrium regime III) at equilibrium: E = 7, μ/E = 1.9, Δλ = 0. (d) Right boundary out-of-equilibrium: E = 7, μ/E = 1.9, Δλ = 1. (e) Left boundary at equilibrium: E = 8, μ/E = 1.7, Δλ = 0. (f) Left boundary out-of-equilibrium: E = 8, μ/E = 1.7, Δλ = 1.

Close modal

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).

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, Tt, and the first assembly time, Tfa (Fig. 5). Specifically, we ran simulations with 105 lattice sweeps (4 × 107 MC steps), with the parameter set {m = 2, 3, …, 9, N = M = 100, N0 = 400, E = 8, μ/E = 1.78, S = 5, Δλ ∈ [0, 1.7]} chosen from the optimal range within regime III.

FIG. 5.

Time spent in the target and first assembly time, normalized by the total number of MC steps. (a) Time spent in the target as a function of the number of target structures. (b) First assembly time as a function of the number of target structures. (c) Time spent in the target as a function of driving force for m = 3. (d) First assembly time as a function of driving force for m = 3. Each marker, “X” or “◦” symbol, represents a result from one realization, and each point represents the median of 14 simulations for (a) and (b) or 18 simulations for (c) and (d).

FIG. 5.

Time spent in the target and first assembly time, normalized by the total number of MC steps. (a) Time spent in the target as a function of the number of target structures. (b) First assembly time as a function of the number of target structures. (c) Time spent in the target as a function of driving force for m = 3. (d) First assembly time as a function of driving force for m = 3. Each marker, “X” or “◦” symbol, represents a result from one realization, and each point represents the median of 14 simulations for (a) and (b) or 18 simulations for (c) and (d).

Close modal

For m > 6, the system does not assemble the desired target structure with or without the driving force, as evident from the values of Tt = 0 [Fig. 5(a)] and Tfa > Ttot [Fig. 5(b)], indicating a storage capacity of mmax = 6 for the parameters chosen.

The trends of Tt and Tfa 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 mmax = 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 Tt) and assembled faster (lower Tfa), 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 Tt [Fig. 5(c)] and the lower values of the first assembly time Tfa [Fig. 5(d)], respectively. However, both Tt and Tfa 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.

Based on nucleation theory, Murugan et al.19 calculated the size of a “critical seed” Sc × Sc for equilibrium self-assembly, below which the system dissolves,

Sc=12μE.
(8)

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 105 lattice sweeps (4 × 107 MC steps), with the parameter set {m = 6, N = M = 100, N0 = 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 Sc = 4.54.

FIG. 6.

Final distance from the desired target structure as a function of seed length. The blue area is the sub-critical seed region. The inset shows the final distance for the maximal equilibrium sub-critical seed, S = 4 < Sc, as a function of the drive value. Each point is the average value over ten simulations, with the error bars representing the statistical error.

FIG. 6.

Final distance from the desired target structure as a function of seed length. The blue area is the sub-critical seed region. The inset shows the final distance for the maximal equilibrium sub-critical seed, S = 4 < Sc, as a function of the drive value. Each point is the average value over ten simulations, with the error bars representing the statistical error.

Close modal

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 pinit = 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 pinit 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 pinit, 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.

Our results show a non-monotonic behavior of the time spent in the target Tt, the first assembly time Tfa, and the final distance df 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 105 lattice sweeps (4 × 107 MC steps), with a similar parameter set {m = 6, N = M = 100, N0 = 400, E = 7, μ/E = 1.78, Δλ ∈ [0, 1.9]}, and examined the characteristic growth time Tg 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.

FIG. 7.

Two quality parameters for varied drive intensity, Δλ: characteristic growth time, normalized by the total number of MC steps (blue) and chimeric value (green). Each “x” or “o” symbol represents one simulation, and each point represents the median of 12 simulations.

FIG. 7.

Two quality parameters for varied drive intensity, Δλ: characteristic growth time, normalized by the total number of MC steps (blue) and chimeric value (green). Each “x” or “o” symbol represents one simulation, and each point represents the median of 12 simulations.

Close modal

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.

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 interact46 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 models20 as the smaller seed can be expanded to the whole target structure by applying the local drive.

Monte Carlo simulations and image interpolation for the heat-maps were run using a Python code, including the libraries NumPy,49 Matplotlib,50 and SciPy.51 

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

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.

The authors have no conflicts to disclose.

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

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

1.
S. A.
Woodson
, “
RNA folding pathways and the self-assembly of ribosomes
,”
Acc. Chem. Res.
44
,
1312
1319
(
2011
).
2.
J. A.
Pelesko
, in
Self Assembly: The Science of Things that Put Themselves Together
, 1st ed. (
Chapman and Hall; CRC
,
2007
), Chap. 3.
3.
K.
Das
,
L.
Gabrielli
, and
L. J.
Prins
, “
Chemically fueled self-assembly in biology and chemistry
,”
Angew. Chem., Int. Ed.
60
,
20120
(
2021
).
4.
E.
Grossman
,
O.
Medalia
, and
M.
Zwerger
, “
Functional architecture of the nuclear pore complex
,”
Annu. Rev. Biophys.
41
,
557
584
(
2012
).
5.
G. C. L.
Wong
,
J. X.
Tang
,
A.
Lin
,
Y.
Li
,
P. A.
Janmey
, and
C. R.
Safinya
, “
Hierarchical self-assembly of F-actin and cationic lipid complexes: Stacked three-layer tubule networks
,”
Science
288
,
2035
2039
(
2000
).
6.
H.
Hess
and
J. L.
Ross
, “
Non-equilibrium assembly of microtubules: From molecules to autonomous chemical robots
,”
Chem. Soc. Rev.
46
,
5570
5587
(
2017
).
7.
P.
Sartori
and
S.
Leibler
, “
Lessons from equilibrium statistical physics regarding the assembly of protein complexes
,”
Proc. Natl. Acad. Sci. U. S. A.
117
,
114
120
(
2020
).
8.
S. F.
Greenbury
,
I. G.
Johnston
,
A. A.
Louis
, and
S. E.
Ahnert
, “
A tractable genotype–phenotype map modelling the self-assembly of protein quaternary structure
,”
J. R. Soc., Interface
11
,
20140249
(
2014
).
9.
Z.
Zhang
and
S. C.
Glotzer
, “
Self-assembly of patchy particles
,”
Nano Lett.
4
,
1407
1413
(
2004
).
10.
A.
Cumberworth
,
A.
Reinhardt
, and
D.
Frenkel
, “
Lattice models and Monte Carlo methods for simulating DNA origami self-assembly
,”
J. Chem. Phys.
149
,
234905
(
2018
).
11.
É.
Fodor
and
M.
Cristina Marchetti
, “
The statistical physics of active matter: From self-catalytic colloids to living cells
,”
Physica A
504
,
106
120
(
2018
).
12.
S.
Ramaswamy
, “
The mechanics and statistics of active matter
,”
Annu. Rev. Condens. Matter Phys.
1
,
323
345
(
2010
).
13.
M. C.
Marchetti
,
J. F.
Joanny
,
S.
Ramaswamy
,
T. B.
Liverpool
,
J.
Prost
,
M.
Rao
, and
R. A.
Simha
, “
Hydrodynamics of soft active matter
,”
Rev. Mod. Phys.
85
,
1143
1189
(
2013
).
14.
C.
Bechinger
,
R.
Di Leonardo
,
H.
Löwen
,
C.
Reichhardt
,
G.
Volpe
, and
G.
Volpe
, “
Active particles in complex and crowded environments
,”
Rev. Mod. Phys.
88
,
045006
(
2016
).
15.
T.
Vicsek
,
A.
Czirók
,
E.
Ben-Jacob
,
I.
Cohen
, and
O.
Shochet
, “
Novel type of phase transition in a system of self-driven particles
,”
Phys. Rev. Lett.
75
,
1226
1229
(
1995
).
16.
E.
Bertin
,
M.
Droz
, and
G.
Grégoire
, “
Boltzmann and hydrodynamic description for self-propelled particles
,”
Phys. Rev. E
74
,
022101
(
2006
).
17.
J.
Madge
and
M. A.
Miller
, “
Design strategies for self-assembly of discrete targets
,”
J. Chem. Phys.
143
,
044905
(
2015
).
18.
J.
Madge
,
D.
Bourne
, and
M. A.
Miller
, “
Controlling fragment competition on pathways to addressable self-assembly
,”
J. Phys. Chem. B
122
,
9815
9825
(
2018
).
19.
A.
Murugan
,
Z.
Zeravcic
,
M. P.
Brenner
, and
S.
Leibler
, “
Multifarious assembly mixtures: Systems allowing retrieval of diverse stored structures
,”
Proc. Natl. Acad. Sci. U. S. A.
112
,
54
59
(
2015
).
20.
W.
Zhong
,
D. J.
Schwab
, and
A.
Murugan
, “
Associative pattern recognition through macro-molecular self-assembly
,”
J. Stat. Phys.
167
,
806
826
(
2017
).
21.
Z.
Zeravcic
,
V. N.
Manoharan
, and
M. P.
Brenner
, “
Colloquium: Toward living matter with colloidal particles
,”
Rev. Mod. Phys.
89
,
031001
(
2017
).
22.
J. D.
Halverson
and
A. V.
Tkachenko
, “
Sequential programmable self-assembly: Role of cooperative interactions
,”
J. Chem. Phys.
144
,
094903
(
2016
).
23.
A.
Murugan
,
J.
Zou
, and
M. P.
Brenner
, “
Undesired usage and the robust self-assembly of heterogeneous structures
,”
Nat. Commun.
6
,
6203
(
2015
).
24.
T.
Mitchison
and
M.
Kirschner
, “
Dynamic instability of microtubule growth
,”
Nature
312
,
237
242
(
1984
).
25.
J.
Boekhoven
,
A. M.
Brizard
,
K. N. K.
Kowlgi
,
G. J. M.
Koper
,
R.
Eelkema
, and
J. H.
van Esch
, “
Dissipative self-assembly of a molecular gelator by using a chemical fuel
,”
Angew. Chem., Int. Ed.
49
,
4825
4828
(
2010
).
26.
M.
Weißenfels
,
J.
Gemen
, and
R.
Klajn
, “
Dissipative self-assembly: Fueling with chemicals versus light
,”
Chem
7
,
23
37
(
2021
).
27.
G.
Makey
,
S.
Galioglu
,
R.
Ghaffari
,
E. D.
Engin
,
G.
Yıldırım
,
Ö.
Yavuz
,
O.
Bektaş
,
Ü. S.
Nizam
,
Ö.
Akbulut
,
Ö.
Şahin
,
K.
Güngör
,
D.
Dede
,
H. V.
Demir
,
F. Ö.
Ilday
, and
S.
Ilday
, “
Universality of dissipative self-assembly from quantum dots to human cells
,”
Nat. Phys.
16
,
795
801
(
2020
).
28.
G.
Bisker
, “
Dissipate your way to self-assembly
,”
Nat. Phys.
16
,
707
708
(
2020
).
29.
A. C. H.
Coughlan
,
I.
Torres-Díaz
,
J.
Zhang
, and
M. A.
Bevan
, “
Non-equilibrium steady-state colloidal assembly dynamics
,”
J. Chem. Phys.
150
,
204902
(
2019
).
30.
C.
Papaseit
,
N.
Pochon
, and
J.
Tabony
, “
Microtubule self-organization is gravity-dependent
,”
Proc. Natl. Acad. Sci. U. S. A.
97
,
8364
8368
(
2000
).
31.
C. B.
Minkenberg
,
L.
Florusse
,
R.
Eelkema
,
G. J. M.
Koper
, and
J. H.
van Esch
, “
Triggered self-assembly of simple dynamic covalent surfactants
,”
J. Am. Chem. Soc.
131
,
11274
11275
(
2009
).
32.
S.
Whitelam
and
R. L.
Jack
, “
The statistical mechanics of dynamic pathways to self-assembly
,”
Annu. Rev. Phys. Chem.
66
,
143
163
(
2015
).
33.
G.
Bisker
and
J. L.
England
, “
Nonequilibrium associative retrieval of multiple stored self-assembly targets
,”
Proc. Natl. Acad. Sci. U. S. A.
115
,
E10531
E10538
(
2018
).
34.
M.
Nguyen
and
S.
Vaikuntanathan
, “
Design principles for nonequilibrium self-assembly
,”
Proc. Natl. Acad. Sci. U. S. A.
113
,
14231
14236
(
2016
).
35.
M.
Tagliazucchi
,
E. A.
Weiss
, and
I.
Szleifer
, “
Dissipative self-assembly of particles interacting through time-oscillatory potentials
,”
Proc. Natl. Acad. Sci. U. S. A.
111
,
9751
9756
(
2014
).
36.
C. J.
Fullerton
and
R. L.
Jack
, “
Optimising self-assembly through time-dependent interactions
,”
J. Chem. Phys.
145
,
244505
(
2016
).
37.
M. E. J.
Newman
and
G.
Barkema
,
Monte Carlo Methods in Statistical Physics
(
Clarendon
,
1999
).
38.
C.-W.
Kok
and
T.
Wing-Shan
, in
Digital Image Interpolation in Matlab
, 1st ed. (
Wiley; IEEE Press
,
2019
), Chap. 2, p.
42
.
39.
M.
Baiesi
and
C.
Maes
, “
Life efficiency does not always increase with the dissipation rate
,”
J. Phys. Commun.
2
,
045017
(
2018
).
40.
E. S.
Andersen
, “
Biomolecular self-assembly: DNA origami rewired
,”
Nat. Nanotechnol.
10
,
733
734
(
2015
).
41.
A.
Tafvizi
,
L. A.
Mirny
, and
A. M.
van Oijen
, “
Dancing on DNA: Kinetic aspects of search processes on DNA
,”
ChemPhysChem
12
,
1481
1489
(
2011
).
42.
M.
Slutsky
and
L. A.
Mirny
, “
Kinetics of protein-DNA interaction: Facilitated target location in sequence-dependent potential
,”
Biophys. J.
87
,
4021
4035
(
2004
).
43.
N. G.
Almarza
and
E. G.
Noya
, “
Phase transitions of a lattice model for patchy particles with tetrahedral symmetry
,”
Mol. Phys.
109
,
65
74
(
2011
).
44.
A.
Reinhardt
and
D.
Frenkel
, “
DNA brick self-assembly with an off-lattice potential
,”
Soft Matter
12
,
6253
6260
(
2016
).
45.
M.
Díez-Minguito
,
J.
Marro
, and
P. L.
Garrido
, “
On the similarities and differences between lattice and off–lattice models of driven fluids
,”
Eur. Phys. J.: Spec. Top.
143
,
269
272
(
2007
).
46.
W. M.
Jacobs
,
A.
Reinhardt
, and
D.
Frenkel
, “
Rational design of self-assembly pathways for complex multicomponent structures
,”
Proc. Natl. Acad. Sci. U. S. A.
112
,
6313
6318
(
2015
).
47.
M.
Sajfutdinow
,
W. M.
Jacobs
,
A.
Reinhardt
,
C.
Schneider
, and
D. M.
Smith
, “
Direct observation and rational design of nucleation behavior in addressable self-assembly
,”
Proc. Natl. Acad. Sci. U. S. A.
115
,
E5877
E5886
(
2018
).
48.
M.
Fialkowski
,
K. J. M.
Bishop
,
R.
Klajn
,
S. K.
Smoukov
,
C. J.
Campbell
, and
B. A.
Grzybowski
, “
Principles and implementations of dissipative (dynamic) self-assembly
,”
J. Phys. Chem. B
110
,
2482
2496
(
2006
).
49.
C. R.
Harris
,
K. J.
Millman
,
S. J.
van der Walt
,
R.
Gommers
,
P.
Virtanen
,
D.
Cournapeau
,
E.
Wieser
,
J.
Taylor
,
S.
Berg
,
N. J.
Smith
,
R.
Kern
,
M.
Picus
,
S.
Hoyer
,
M. H.
van Kerkwijk
,
M.
Brett
,
A.
Haldane
,
J. F.
del Río
,
M.
Wiebe
,
P.
Peterson
,
P.
Gérard-Marchant
,
K.
Sheppard
,
T.
Reddy
,
W.
Weckesser
,
H.
Abbasi
,
C.
Gohlke
, and
T. E.
Oliphant
, “
Array programming with NumPy
,”
Nature
585
,
357
362
(
2020
).
50.
J. D.
Hunter
, “
Matplotlib: A 2D graphics environment
,”
Comput. Sci. Eng.
9
,
90
95
(
2007
).
51.
P.
Virtanen
,
R.
Gommers
,
T. E.
Oliphant
,
M.
Haberland
,
T.
Reddy
,
D.
Cournapeau
,
E.
Burovski
,
P.
Peterson
,
W.
Weckesser
,
J.
Bright
,
S. J.
van der Walt
,
M.
Brett
,
J.
Wilson
,
K. J.
Millman
,
N.
Mayorov
,
A. R. J.
Nelson
,
E.
Jones
,
R.
Kern
,
E.
Larson
,
C. J.
Carey
,
İ.
Polat
,
Y.
Feng
,
E. W.
Moore
,
J.
VanderPlas
,
D.
Laxalde
,
J.
Perktold
,
R.
Cimrman
,
I.
Henriksen
,
E. A.
Quintero
,
C. R.
Harris
,
A. M.
Archibald
,
A. H.
Ribeiro
,
F.
Pedregosa
,
P.
van Mulbregt
 et al, “
SciPy 1.0: Fundamental algorithms for scientific computing in python
,”
Nat. Methods
17
,
261
272
(
2020
).

Supplementary Material