Active motion of synthetic nanomotors in ﬁlament networks

The collective behavior of chemically powered nanomotors in complex ﬁlament networks is investigated. The synthetic motors are small oligomers and use self-diffusiophoresis for their active motion. Much like biological molecular machines, these motors attach to the ﬁlaments, which highly suppress orientational Brownian motion, and move along them. The collective motion is inﬂuenced by structure of the network and the forces that bind the motors to the ﬁlaments; it is characterized by strong clustering, especially near cross-links in the network, and differs substantially from that in a simple ﬂuid phase where only weak clustering occurs. The features of the collective behavior are quantitatively examined through computations of the effective binding potentials and the dynamics of probe oligomers. The results are relevant for applications in the biological and materials sciences involving complex natural and synthetic networks activated by small chemically powered motors.


I. INTRODUCTION
Molecular motors with linear dimensions of tens of nanometers carry out their functions in complicated and crowded environments.In the cell, motor proteins, such as those in the family of kinesins, actively transport material along microtubules in the cellular medium comprising fluid and a variety of small molecules and large macromolecules within a dense network of filaments of different types [1][2][3][4].These and the many other kinds of molecular motors derive their fuel from an equally complex network of biochemical reactions, the whole operating under nonequilibrium conditions.Despite this, complexity biological systems function effectively, even in the presence of strong thermal fluctuations.On larger micrometer scales, bacteria and other microorganisms are able to swim in environments, such as gel media, and in microfluidic structures [5,6].
This paper deals with the behavior of synthetic motors that also use chemical energy as their fuel source in systems with filament networks that mimic those in the cell or gel media.An understanding of how these motors function under such complex environmental conditions is required to exploit their full potential in diverse biological and materials science applications.In contrast to molecular machines that typically use conformational changes to effect directed motion, the motors in our paper do not rely on moving parts for their propulsion and, instead, operate by a self-diffusiophoretic mechanism [7][8][9][10].In diffusiophoresis, a colloid in a chemical concentration gradient induces fluid flows in its vicinity that lead to motion in media free from external forces [11][12][13].
Extensive research on these objects has resulted in the construction of numerous different motors with sizes ranging from tens of nanometers to several micrometers using various types of chemical fuel and operating by different mechanisms [14][15][16][17][18][19][20][21].The behavior of large collections of such motors has been frequently a topic of current research [22][23][24][25][26][27][28][29][30][31][32][33].Although many of these investigations consider motor dynamics in simple fluid-phase systems or in the presence of walls [34], there have been investigations of motors in more complex media including those with inert and active obstacles as well as chemical patterns [35][36][37].
Molecular motors mediate some of the effects of thermal fluctuations and carry out directed transport by attaching to the macromolecular filament network in the cell.Similarly, if the medium in which the synthetic chemically powered motors operate contains macromolecular filaments, the motors may attach to the filaments and move along them [38].This attachment serves to "tame" the orientational Brownian motion that affects such motors in solution, which is especially important for small nanoscale motors where reorientation is rapid.When there are many motors in the system and the filament network has a complicated structure, we show below how the collective motor dynamics is modified.
The outline of the paper is as follows.Section II describes the motor used in this paper, how the filaments and filament network are constructed, and the dynamics used to evolve the C N N A B FIG. 1.A chemically powered oligomeric motor.The motor consists of three beads: one catalytic bead (C) and two noncatalytic beads (N).A catalytic reaction occurs when fuel particles (A) encounter the C bead and are converted to product particles (B).
entire system.In Sec.III, in order to extract some general features of the dynamics of motors on filaments in a simple context, the properties of the diffusiophoretic motors in bulk solution is contrasted with those when there is just a single filament or two cross-linked filaments in the system.The full complexity of the behavior of motors in filament networks is the subject of Sec.IV where motor and filament network properties are analyzed in detail.The conclusions are given in Sec.V.

II. MOTORS, FILAMENT NETWORK, AND DYNAMICAL MODEL A. Motors, filaments, and filament networks
The synthetic self-propelled motor we consider is a stiff short chain made of three linked beads; one end bead is catalytic (C) whereas the remaining two beads are noncatalytic (N) (Fig. 1).The three-bead oligomeric motors can attach to a filament and move along it whereas their orientational motion is highly suppressed [38].All motor beads have the same radii σ C = σ N = 2.0.(All results and parameters are given in the dimensionless units defined in the Appendix.)The beads are connected by stiff harmonic springs V bond = 1 2 k s (r − r eq ) 2 where the spring constant is k s = 500 and the equilibrium bond length is r eq = 3.0, whereas the bending stiffness is taken into account by adding a spring linking the end beads with 2r eq .The volume of the motor is computed by subtracting the overlap volumes between the beads, and the total mass of the motor is m M = 848.23.
The filaments are linear chains made from 50 beads (F), each with radius σ F = 0.75.The neighboring beads in a chain are linked by stiff harmonic springs with the same forms as those for the motor but with an equilibrium bond length of r F eq = 1.0, so the length of a filament is F = 50.The bending stiffness of a filament is controlled by a three-body potential, V bend = k b [1 − cos θ ] with k b = 100, and cos θ = ri−1,i • ri,i+1 , where ri, j = (r i − r j )/|r i − r j |.The total mass of the filament is m F = 795.21,accounting for overlap of the filament beads.
A complex semiflexible network is built from N F = 50 cross-linked filaments of the same length F as follows.The N F filaments are first randomly distributed in space to ensure that the network is homogeneous, and then cross-links are added to provide structural stability to the network.The characteristics of the network are determined by the mesh size of the network ξ or the average distance between crosslinks c .The filaments cannot form loops or knots since their persistence length is much longer than the characteristic mesh size.Some beads in each filament are chosen to be cross-linker beads with the average distance between cross-links given by c .The cross-links are stable because the potential between cross-linker beads α and α in different neighboring filaments is attractive and is given by [39] V att (r where FF = 10, σ FF = 2 σ F + σ with σ = 1, r c = 2 1/6 σ FF , and w αα = 2σ F .A example of a filament network constructed in this way is shown in Fig. 2.

B. System dynamics
The motors and filaments are contained in a cubic simulation box of linear size L = 50 with periodic boundary conditions in all directions.The fluid phase contains N = N A + N B particles of species A and B with mass m = 1.With the exception of the harmonic spring potentials discussed above and used to construct the motors and filaments, all other intermolecular interactions take place through repulsive Lennard-Jones (LJ) potentials of the form where θ (r) is the Heaviside function and the separation between a particle i of type α and a particle j of type α is Interactions among fluid particles are described by multiparticle collision dynamics [40][41][42], and the time evolution of the entire system is carried out using a method that combines molecular dynamics with multiparticle collision dynamics for solvent particles.Our simulations are carried out under low Reynolds number conditions where inertia plays no role.Additional details of the simulation method are given in the Appendix.
In propulsion by self-diffusiophoresis chemical reactions on the catalytic portions of the motor lead to an inhomogeneous distribution of reactants and products in the vicinity of the motor.If these inhomogeneously distributed fluid species interact with the motor through different intermolecular potentials, they give rise to forces on the motor which, by momentum conservation, lead to fluid flows in the surrounding medium that are responsible for propulsion [8].In the oligomeric motor considered here, the reversible chemical reactions A + C B + C that power motor motion take place on the catalytic bead of a motor.Since self-propulsion cannot take place at equilibrium, for the mechanism to operate, the system must be maintained in a nonequilibrium state where detailed balance is broken.
The reversible reactions microscopically in our simulations through collision events between the A or B particles and the C motor beads where a transformation from one species to the other take place with unit probabilities p + = p − = 1 when the reaction radius R c is crossed by an A or B particle [43].The radius R c is taken to be infinitesimally greater than the range of the potentials so that there are no abrupt changes in the forces.These collision events determine the reaction rate coefficients k ± whose values depend on the degree to which the rates are limited by reaction on the catalytic surface or diffusion to it.Assuming that the entire surface of the catalytic sphere is available for catalytic reactions, approximate values for these rate coefficients can be obtained from kinetic theory which yields [44] , the Smoluchowski diffusion controlled rate coefficient with D M and D, respectively, the motor and common A and B species diffusion coefficients.From simple collision theory, the intrinsic rate coefficient is k is the reduced mass.In the simulations, the motor mass is much larger than the reactive species mass m M m, also D M D. For our parameters, the reaction rate is strongly influenced by diffusion but is neither fully diffusion nor reaction limited.Nonreactive collisions of the A and B species with the noncatalytic beads are governed by the intermolecular interactions described above.
The reactive collision events satisfy detailed balance and, starting from an initial state, the system as described above will evolve to equilibrium where k + /k − = c eq B /c eq A = K eq with c eq α as the equilibrium concentrations of the species α = A, B, and active motion will cease.( Since we have chosen p ± = 1, k + = k − , and the equilibrium constant K eq for this reaction is unity.) The system can be maintained in a nonequilibrium state by coupling to reservoirs with fixed concentrations of chemical species or, as in biochemical systems, by reactions in the environment that are themselves maintained in nonequilibrium conditions.Here, we impose nonequilibrium conditions through fluid phase reactions B k +b k −b A which break detailed balance [43] and occur locally in the bulk of the solution with rate coefficients k ±b .These fluid phase reactions were implemented using the reactive version of multiparticle collision dynamics [45].Detailed balance is broken because the fluid phase reaction occurs under nonequilibrium conditions by a different mechanism where the constant concentrations of pool chemical species (not shown) are incorporated in the k ±b rate coefficients [43].The deviation from equilibrium is controlled by the (dimensionless) chemical affinity [10], where detailed balance will be broken, and active motion will persist.
The propulsion velocity of the oligomeric motor is directed along its long axis, V sd = V sd ûM , where ûM is a unit vector directed from the noncatalytic to catalytic beads and V sd and can be written in the general form [10] where W rxn is the reaction rate, which is proportional to the chemical affinity A rxn and χ is the diffusiophoretic parameter that depends on the fluid and motor properties including the interaction potentials.We introduce the following terminology: A forwardmoving motor moves with the catalytic bead at its head whereas a backward-moving motor moves with the last noncatalytic bead at its head.On the basis of Eq. ( 4), the motor can be made to move forward or backward by changing either the chemical affinity or the interaction potentials.For a given set of interaction potentials, changing the ratio k +b /k −b from greater to less than unity changes the sign of the affinity and the direction of motor motion.
Differences in the potential energies can be due to interactions of any or all fluid reagent species with either the catalytic or noncatalytic portions of the motor.We have chosen the fluid-motor interaction potentials such that the potential functions for the A and B species differ only in their interactions with the noncatalytic motor beads.In this circumstance, the catalytic bead is responsible for producing the concentration inhomogeneities, and forces due to the noncatalytic beads lead to propulsion.For the choice in this paper, ( NA = 1.0,NB = 0.1) and positive chemical affinity (k +b > k −b ), the motor moves in the forward direction due to the stronger repulsive forces at the noncatalytic end of the motor.If these interaction energies were reversed ( NB = 1.0,NA = 0.1), the motor would move backward for the same value of the affinity.A description of the dynamics for a wide range of possible interaction potential differences has been given for Janus colloids [40].In this paper, we fix the interaction potentials and control the direction of motor motion by the chemical affinity.
where N F and N M are the numbers of filaments and motors, respectively.The total simulation time is denoted by t m (∼ 10 5 ).

A. Motor dynamics in the fluid phase
In order to put our study of motor dynamics in filament networks into perspective, we first briefly describe the dynamics of the three-bead oligomeric motors in solution.Under nonequilibrium conditions, the motor will execute directed motion along its symmetry axis ûM , a unit vector directed from the noncatalytic to catalytic beads in a motor.
As discussed above, the nonequilibrium conditions in the system are specified by the values of the fluid-phase reaction rate coefficients.Three different choices of these fluidphase rate coefficients are considered: (i) (k +b = 10 −2 , k −b = 10 −3 ) where the motor moves in a forward direction (f-m); (ii) (k +b = 10 −3 , k −b = 10 −2 ), the motor moves in a backward direction (b-m); and (iii) (k +b = k −b = 10 −3 ) where the system satisfies detailed balance and leads to an equilibrium state (inactive).Note that the directed motion of the motor, forward or backward, is not controlled by changes in the interaction potentials of the reactive species with the motor beads; these remain unchanged in our paper.
The self-diffusiophoretic propulsion velocity of the center of mass of the motor V sd can be determined in the simulations by computing where V is the velocity of the center of mass of the motor and • • • denotes the average over time and realizations.
In addition to directed motion along ûM , an active motor will undergo orientational Brownian motion that leads to longtime diffusive dynamics characterized by an effective diffusion coefficient where the orientational relaxation time τ R can be determined from the decay of the orientational correlation function, ûM (t ) • ûM (0) = e −t/τ R .
Table I(a) shows how the average motor propulsion velocity of single active motor V sd , the orientational relaxation time τ R , and the effective diffusion coefficient change with the values of k +b and k −b .For reference, a single inactive motor in an equilibrium system with (k +b = k −b = 10 −3 ) has dif- fusion coefficient D M = 0.0035 and orientational relaxation time τ R = 1400.Diffusion is anisotropic for the three-bead oligomeric motor, and D M is the average of its parallel and perpendicular components Next, we consider the collective dynamics of these oligomeric motors in solution.Depending on the nonequilibrium conditions, dynamic clusters with a finite lifetime can occur.The tendency to form small transient clusters can be quantified by considering the steady-state radial distribution functions g(r) with r being the magnitude of the distance between the center spheres of the motors, which is defined as where V is the volume of the system, N M is the number of motors, and the angle brackets again denote an average over time and realizations.Figure 3(a) shows plots of this radial distribution function.From the presence of a peak in g(r), one can see that small transient clusters are formed for forward-moving motors whereas there is no such tendency for backward-moving and inactive motors.
The greater tendency of forward-moving motors to form clusters is similar to that for self-diffusiophoretic motors in fluids.For our parameters, forward-moving motors move towards higher product B concentrations, whereas backwardmoving motors move towards lower B concentrations.Consequently, in many-motor systems, a forward-moving motor will detect a higher-B concentration from a nearby motor and move towards it, but a backward-moving motor will tend to avoid other motors because of the B-concentration field they produce.
Additional insight into this tendency to self-assemble can be gleaned from an examination of the local spatial orientational order described by the function, (7) where is the number density of motor pairs separated by the distance r.From this definition, two motors separated by the distance r have the same or opposite orientations when C u (r) > 0 or C u (r) < 0, respectively.Correspondingly, in Fig. 3(b), one sees that the main distinctive features of the curves are in their structures at small separations r < 7.5.Since C u (r) < 0 in this range of r values, nearest-neighbor forward-moving motors tend to have an antiparallel alignment, whereas neighboring backward-moving motors have C u (r) > 0, and they tend to have a parallel alignment.
The orientational trends have their origin in the clustering tendencies discussed above.In our motors, it is the noncatalytic beads that respond to the concentration gradients; thus, for forward-moving motors noncatalytic beads of one motor will tend to be attracted to the catalytic bead of the other motor so that an antiparallel configuration is favored.The opposite is true for backward-moving motors.Since their noncatalytic beads tend to avoid the catalytic beads, a parallel arrangement is favored.
Since detailed balance is not violated when k +b = k −b , the chemical reactions on the catalytic spheres of inactive motors take place at equilibrium.The interaction potentials remain the same as for the active motors, and the figure shows the much smaller equilibrium spatial and orientational structural features of a suspension of inactive oligomers in a fluid with these interactions.

B. Motor dynamics when attached to filaments
The dynamics described above changes if a filament is present in the system.Although the motor and filament beads interact directly through repulsive interactions, solventmediated effects give rise to effective attractive interactions that lead to binding of motors to filaments.The strength of these effective attractive interactions will determine length of time motors remain attached, and the corresponding average fractions of bound and unbound motors in the system under steady-state conditions.For strong binding, motors will tend to move along filament for long periods of time, and the orientational Brownian motion will be suppressed leading to more efficient directed motion.
The changes in motor properties that occur as a consequence of attachment to a filament can be seen in a comparison of results in parts (a) and (b) in Table I.The magnitudes of motor velocities are somewhat larger but within statistical uncertainties, and the motor reorientation times are substantially longer when motors are attached to a filament.
To investigate and distinguish the single motor behavior on different timescales, the mean-squared displacement (MSD) r(t 2 of the motor is shown in Fig. 4. In bulk solution, the active and inactive motors display the expected behavior: The active motor has a distinct ballistic regime followed by a long-time diffusive regime, whereas the inactive motor has a very short inertial regime followed by diffusive dynamics.By contrast, for active motors on the filament, ballistic motion persists for long times, and the diffusive regime is not reached on the timescale of the simulations τ R > t m ∼ 10 5 .The diffusive dynamics on the filament for the inactive motors is characterized by a larger value of the diffusion coefficient D 0 = 0.0063, compared with that for inactive motors in the solution D M = 0.0035 as mentioned in Sec.III A. FIG. 4. Log-log plots of the mean square displacement r 2 (t ) for active and inactive motors in bulk solution and attached to a single filament.Fits to the mean-square-displacement curves are indicated by dashed lines for the ballistic ∼V 2 u t 2 and dotted lines for the diffusive ∼t portions of these curves.Only the ballistic regime is seen for motors on a filament on the timescales of the plot.
Since the filament networks we will consider are crosslinked, it is useful to examine how motors interact with the cross-linked segments of the filament network.For this purpose, it is sufficient to consider two pinned cross-linked filaments which, for simplicity, we take to be linked at right angles to each other.The presence of such a crosslink can substantially alter the effective attractive interactions between the motor and the filaments.This effect can studied quantitatively through the effective potential of mean force that the motor experiences in the vicinity of the filaments, which is generally defined in terms of the radial distribution function by We first computed the potential of mean force W 1 (r) = −k B T ln g MF (r) for the interaction of the central motor bead with a bead in a single filament where no cross-links are present.Here, the radial distribution function g MF (r) is defined as in Eq. ( 6) with r, now, the distance between the central bead of a motor and a filament bead.Since an active motor can move along the filament, we used umbrella sampling [46] with a biasing harmonic potential to localize the motor to a specific bead in the filament.Figure 5(a) is a plot of W 1 (r) versus r and shows an attractive interaction with a well depth of approximately 0.6 (3.0k B T ).For two perpendicular cross-linked filaments, we determined the effective potential W 2 (r) = −k B T ln g MC (r) between the central bead of the motor and the midpoint between the two cross-linker beads.Once again, the radial distribution function g MC (r) was determined using umbrella sampling with a harmonic biasing potential that localizes the motor to the vicinity of the cross-linker midpoint.The results of this calculation are shown in Fig. 5(b).Now, one sees a doublewell attractive potential with well depths of approximately 1.3 (6.5k B T ) and 1.0 (5.0k B T ), separated by a high barrier.The double-well structure of this effective potential has its origin in the two predominant orientations the motor can take when in the vicinity of the cross-linked filaments: (1) The motor is attached to and lies parallel to one of the filaments so that its center is in close proximity to the filament linker beads [left inset of Fig. 5(b) corresponding to the short-distance minimum], and (2) the head of an attached motor encounters the other filament at their intersection point and, in this case, the distance between the center bead of the motor and the midpoint of the linker beads is r ≈ 5.6 [right inset Fig. 5(b) corresponding to the longer-distance minimum).The motor can switch between these configurations, but such switches are rare events since the free-energy barrier is approximately 5 to 6k B T .
The dynamics changes for systems with two cross-linked filaments and many (20 are considered here) motors.After a short transient time period, all 20 motors attach to and remain on the filaments for the duration of the simulation.Motors are observed to move along the filaments, but those attached to the cross-links remain there for long times and hinder the movement of other motors, leading to the formation of dynamic clusters around the cross-linked region.It follows from this that the magnitude of average motor velocity is much smaller than that for a single motor as can be seen in Table I The enhanced binding of motors to cross-links will influence the collective dynamics when motors are attached to filaments.To examine this effect in a simple context, we present results for the potential of mean force W 2 (r) for a system with 20 motors and two orthogonal cross-linked filaments in Fig. 6. Results for motors moving forward and backward as well as inactive motors with no net propulsion are shown in the figure.The double-well structure persists but is altered, and additional structure at longer distances appears as a result of clustering near the cross-linkers.Note that the most prominent double-well structure is seen in the plot for inactive motors.This is likely due to their preferential attachment to the cross-linkers and their lack of a diffusiophoretic force to assist them to overcome such potential traps.Since backwardmoving motors do not favor cluster formation, they are able to escape the trapping regions more effectively than inactive motors, and this is reflected in the smaller depths of the double wells and the added structure at longer distances.Forward-moving motors show the most complicated structural features in W 2 (r).The potential well at short distances has a depth comparable to that for inactive motors, but a much more complex structure is seen at longer distances.Since forward-moving motors tend to form transient clusters, and these clusters lie predominantly near the linker points, this favors the formation of a strong short-range attractive potential.However, because the motors are active, they are able to explore a larger configuration space, which leads to the additional structural features seen in the plot.
We may also compare these forward-moving motor results for a system with 20 motors with those discussed earlier for a single motor in Fig. 5(b).In the many-motor system, the simple double-well structure is highly modified, and the well depths are much smaller since many motors accumulate in the limited space around the cross-link.

IV. COLLECTIVE MOTOR DYNAMICS IN FILAMENT NETWORKS
The semiflexible filament networks are constructed using the method given in Sec.II and as discussed there for networks with cross-links, the mesh size ξ or average distance between cross-links c depends on the lengths of the filaments F and the distribution of the cross-links.The simulations in this section involve networks constructed from N F = 50 filaments each with length F = 50, and N c ∼ 90-110 crosslinks separated by an average distance of c ∼ 4-8.The initial state of the system is constructed from a realization of the filament network, fluid particles, and N M = 45 chemically active motors that are randomly placed in the fluid phase.

A. Macroscopic filament properties
Before describing the collective dynamics of motors, it is first useful to quantitatively assess the physical characteristics of the filaments in a network with cross-links with and without the presence of motors.For this purpose, the filament bending modulus was determined from an analysis of the thermal fluctuations of the filament chain conformations.A model for the force-extension response of such semiflexible filaments is based on an energy functional that accounts for the bending energy and tension [39,47].It takes the form Here, u(z) is the transverse deviation of the filament away from the line joining the end points of the filament with contour length F , u is the local curvature, and u is the local axial strain along the filament.Letting ni , i = 1, 2 denote unit vectors in the two transverse directions, we can write u(z) = [u 1 (z), u 2 (z)].The components u i (z) can be expressed in a Fourier series as u i (z) = q u i (q) sin(qz) with wave-numbers q = π/L, 2π/L, . ... The Fourier coefficients u i (q) in the two transverse directions satisfy the equipartition theorem; therefore, the resulting power spectrum of these fluctuations in either transverse direction is given by [39,47] where κ is the filament bending modulus and τ is its tension.Figure 7 shows simulation results for the power spectra of systems of filaments with and without cross-links as well as filament networks with cross-links and active motors and compares them to S(q) obtained using Eq. ( 10).In the simulations, S(q) was obtained by fast Fourier transform of the filament fluctuation amplitudes in the two transverse directions and averaging over all filament configurations.
From an examination of this figure, one sees that the power spectrum scales as S(q) ∝ q −4 for small values of the wave-number q.The bending rigidities were obtained FIG. 8.A instantaneous configuration of chemically powered motors on the cross-linked filament network.by fitting Eq. ( 10) to the simulation data in the bendingdominated regime.For a filament whose length F is less than the persistence length p , F p = κ/k B T , the filament is nearly straight with only small transverse fluctuations.With no cross-links, the filaments are stiff since they were designed to have a bending parameter of k b = 100, and the power-law fits of the data to Eq. ( 10) do, indeed, yield κ ∼ 100 and γ = 0.This result indicates that the bending stiffness of a filament network with no cross-links is the same as that of an individual filament.
The presence of cross-links gives rise to strong interactions among filaments which no longer behave as individual elements but form part of an interconnected network with different properties.For such a cross-linked network, we find that the average bending rigidity of the filaments in the network decreases to κ = 50 since the forces from the other filaments tend to reduce the filament stiffness.If, in addition, motors are attached to the filaments in the network the bending rigidity decreases further to κ = 40.This smaller change indicates that the bending rigidity, which is an essential feature of the filament network, is mainly governed by the cross-linked structure of the network.

B. Collective dynamics in the filament network
Building on the two-filament results presented in Sec.III, we now consider the dynamics of motors in a complex semiflexible filament network.Starting from an initial state where the motors are placed uniformly at random in the system, and accounting for exclusion by the filaments, the motors quickly attach to the filaments.An example of a configuration showing motors on a cross-linked filament network formed by such a process is given in Fig. 8.
The collective behavior of motors in systems with a high concentration of motors per unit filament length differs substantially from that in the bulk phase or in networks with low motor densities.Table I (d) shows that the motor propulsion velocity decreases by more than an order of magnitude for this (N M = 45, N F = 50) system, compared to the results in (a)-(c) in this table.
The position dependence of the effective interaction W eff between the central motor beads and the midpoints between the cross-linked filament beads can again be used to characterize the clustering properties of motors on filaments.Since the links are uniformly distributed and the motors form clusters around the links, it is difficult for motors to escape the trapping regions around the links.The effective potentials W eff (r), analogous to those in Fig. 6, are shown in Fig. 9.One sees that the potential curves for forward-moving, backward-moving, and inactive motors are almost indistinguishable.The two dominant minima are in approximately the same locations as in Figs.5(b) and 6 but are shifted slightly to shorter distances; however, the well depths are significantly smaller.The average over many cross-links in the network tends to eliminate the finer structure seen in Fig. 6 where a single cross-link is considered so that only the structure corresponding to the dominant structures survives.From these results, one can conclude that the steady-state structural features of the motor clusters are governed mainly by the structure of the filament network and especially the density of the cross-links.

C. Filament network dynamics
The static structural properties of the motor-filament network were examined in Sec.IV A. Here, we investigate whether active motors can change the dynamical properties of the filament network.In this connection, we study a specific aspect of the dynamics: The changes in the forces exerted on probe particles by the network.The probe particles we consider are three-bead inactive oligomers which are tethered to randomly chosen locations in the system by harmonic springs where k p is the spring constant and r c0 is the equilibrium separation.The tethering points do not lie on the filaments but, such as the untethered motors, the tethered oligomers can attach to the filaments and, thus, serve as probes of the network dynamics by monitoring the forces they experience through such probe-filament interactions.The effects of these interactions are reflected in the MSD of the tethered oligomers, which behave much like probe molecules confined by optical traps.
In the absence of any interactions with the filaments or motors, a probe oligomer's overdamped dynamics is well described by the Langevin equation, where the random force satisfies the fluctuation-dissipation relation f p (t ) f p (t ) = 2D p 1δ(t − t ) with D p = k B T /ζ p as the diffusion coefficient and ζ p is the friction coefficient.The MSD of the probe oligomer computed from this equation is This expression accurately describes the time evolution of the simulation data with a value of ζ p = 60 for the probe oligomer friction coefficient as can be seen in the inset of Fig. 10.As expected, the MSD saturates at a value of 6k B T /k p = 4.8 when k p = 0.25 and at half the value if k p = 0.5.As mentioned above, the probe oligomers attach to the filaments, and the effects of forces exerted on the tethering springs can be monitored through the mean-squared displace-TABLE II.Dynamical properties of the probe oligomers for various equilibrium and nonequilibrium conditions.The effective diffusion coefficients Dp , the effective spring constant kp , and the friction coefficient ζp from fits of the simulation data to Eq. ( 13  ment r(t ) 2 .Provided the number of such probe oligomers is small, they will not modify the properties of the network.Below, we consider N p = 5 probe oligomers; this number was tested and found to satisfy this criterion.The sensitivity of these probes can be tuned by changing their force constants k p .
The results of calculations of r(t ) 2 are shown in Fig. 10 for the filament network in the absence of motors and with forward-moving motors, backward-moving motors, and inactive motors.The temporal behavior of r(t ) 2 in the presence of the filament network differs qualitatively from that for a probe oligomer in solution.The most obvious qualitative change in panel (a) of this figure for k p = 0.5 is the continual linear increase in r(t ) 2 up to the times t ≈ 10 4 -10 5 considered in our simulations.In panel (b), where the spring constant is smaller k p = 0.25, the linear growth regime is now resolved, and one can see the differences in behavior when motors are active in the filament network.
The phenomenological expression, can capture the gross structure of the probe oligomer MSD curves.It shows that the short-time behavior is similar to that of a probe oligomer in solution but with modified force constant and friction to account for the interactions with the network and motors.The friction coefficient ζp is much larger than its solution value and gives rise to the longer times needed to reach the diffusive regime and a plateau value.In this intermediate-time region, the MSD curves for the active forward-and backward-moving motors are higher than those with inactive and no motors due to changes in the effective force constants as seen in Table II.The values of ζp and kp were determined from fits to the short-time regime of the MSD curves in Fig. 10, whereas Dp was computed from the slope of a linear fit to the data in the longer-time linear regime.Since the probes are strongly bound to the massive filament and motor network, the linear diffusive regime characterized by small effective diffusion coefficients Dp (see Table II) reflects the diffusive dynamics of the entire motor-filament network and its internal dynamics.Due to the harmonic nature of the tethering springs, r(t ) 2 must saturate at sufficiently long times, but our simulations have not yet reached this regime.The damped oscillatory behavior seen in some of the curves is likely due to the oscillator regime change from strongly overdamped dynamics to more weakly damped dynamics because of coupling to the network.
From these results, we see that the gross changes to the probe MSDs, in the presence of the filament network, are the effective force constants and friction coefficients that describe the short-time dynamics and the long linear increase in the MSD curves before saturation at much longer times (not seen on the timescales of the simulation).These major changes depend primarily on the transient binding of the probes whether active or inactive to the filaments as evident in the data for the larger force constant in Fig. 10(a) where the MSD curves for forward, backward, and inactive motors are very similar.Effects that are specific to the different motor types are more nuanced and are seen only if the probe particles are tethered by weaker harmonic springs [Fig.10(b)].Thus, the network structure strongly influences the motions of the diffusiophoretic motors in this paper, but motor activity only weakly changes the network dynamics.However, if the filaments are less rigid, the motor motion may couple more strongly to internal network collective modes.

V. CONCLUSION
In solution, small oligomeric motors execute directed motion along their long axis and are subject to strong thermal fluctuations that lead to rapid orientational motion.This relatively rapid thermally driven orientational motion restricts the timescale on which the active ballistic directed motion occurs, and strongly influences the character of the collective dynamics in many-motor systems.
By contrast, in complex filament networks, these motors have a propensity to attach to the filaments and move along them, much like molecular motors in the cell.Orientational Brownian motion is highly suppressed giving rise to longlived active ballistic motion.Cross-links in the filament network play important roles in determining the characteristics of the collective dynamics.The combination of stronger binding forces in these regions and the fact that crossing filaments can act as obstacles to motor motion lead to strong clustering in these regions of the network.The effects are very strong for forward-moving motors that have stronger clustering tendencies even in the bulk phase.The use of probe oligomers to investigate aspects of the filament network dynamics shows how the character of motor motion, forward or backward, is reflected in form of the mean-square displacement of the probe oligomers.
The scales on which the phenomena described in this paper can be observed are diverse.Thick biofilaments, such as microtubules have diameters in the range of 25 nm.Given the relative sizes of the motors and filaments in our paper, this corresponds to motors with diameters of approximately 70 nm with lengths of about 200 nm, comparable to those of some of the very small motors that have been constructed and studied experimentally [19,48].This indicates that biological applications may be possible for synthetic motor-biofilament systems.On larger scales, synthetic filament networks with much thicker filaments have been constructed using threedimensional printing methods [49], and such networks activated by micrometer-sized chemically powered motors have the potential for materials science applications.
The work described in this paper indicates that the ability of small synthetic motors to perform useful transport tasks in filament networks depends on a number of factors; some are similar to those for molecular machines whereas others are unique to the synthetic motors considered here.In qualitative terms, molecular machines, such as the kinesins walk along filaments in specific directions determined by the filament polarity and when not bound undergo diffusion in the medium.They typically remain bound to filaments for times that are sufficiently long so that many steps along the filament can be taken for efficient cargo transport.The diffusiophoretic motors considered here propel themselves in bulk solution along their polar axis determined by the location of the catalytic region.If the motors are small, they will reorient rapidly leading to enhanced diffusion but with no specific target direction.We have seen that attachment to filaments can tame this orientational motion, but if the motors are strongly bound to the filaments (this binding is especially strong near cross-links), their propulsion velocity is smaller; thus, the resulting dynamics depends on competing effects.Since the parameter space that fully characterizes the filament network, motor dynamics, and binding properties is large, these parameters can be used to design motor-filament systems that have desirable properties.For example, if the binding of motors to the filaments is much weaker, the active motors may be able to surmount free-energy barriers that lead to strong clustering at cross-links.In this way, motor-filament interactions may be used to control the average length of time the synthetic motors remain attached to filaments.Then, such as molecular motors, they will attach and detach in the course of their motion on a filament.In a filament network, a motor that detaches from one filament may attach to another filament and continue its active motion.The motors are still able to form clusters in regions with cross-links.
These are not the only aspects of systems that merit study.Our motors move equally well in either direction along filaments, but it is possible to design motors with biased motion.In addition, it has been observed in our studies that motors can actively transport small filaments through the network, much like the active transport of organelles and other materials by molecular machines in the cell.These and other aspects make synthetic motor motion in complex networks an interesting topic to study.

FIG. 2 .
FIG. 2. A filament network made from N F = 50 cross-linked filaments.The beads of the filament chains are rendered in cyan, whereas the cross-linker beads are yellow.
r i j = |r i − r j |.The interactions between fluid particles and motor and filament beads have LJ energy parameter αS and cutoff distance r c = 2 1/6 σ αS , where σ αS = (σ α + σ S ).The repulsive potential between two beads in different motors has σ MM = 2σ C + σ with σ = 1.Interaction strengths among motor and filament beads are MF = 0.1, and among motor monomers, MM = 1.0.Interaction strengths between solvent particles and filament beads are AF = BF = 0.1.Filament beads interact with those in neighboring filaments with strength FF = 1.0.The A and B fluid particles have identical effective radii σ A = σ B = 0.25 and energy parameters AC = BC = NA = 1.0,NB = 0.1 for their interactions with the motor beads.

FIG. 5 .
FIG. 5. (a) The potential of mean force W 1 (r) between the single forward-moving motor and the filament N M = 1, N F = 1.(b) The potential of mean force W 2 (r) between a single forward-moving motor and the midpoint between two cross-linker beads in a system with two perpendicular cross-linked filaments N M = 1, N F = 2.The results were obtained from time averages over 20 independent realizations of the dynamics.

FIG. 7 .
FIG. 7. The power spectrum of the thermal filament fluctuations in the filament network with N F = 50.The bending rigidity of a single filament is determined by k b = 100.The bending rigidities are obtained from fits of Eq. (10) (solid lines) to the simulation data (solid circles) as shown in the figure: (i) κ = 100 for the filament network without the cross-links (red solid circles); (ii) κ = 50 for the filament network with the cross-links (green solid circles); (iii) κ = 40 for the cross-linked filament network with N M = 45 active forward-moving motors (blue solid circles).

TABLE I .
Dynamical properties of a motor for various equilibrium and nonequilibrium conditions: The table gives the average motor propulsion velocity V sd , the orientational relaxation time τ R , and the effective diffusion constant D eff for different numbers of motors N M and different numbers of filaments N F