# Unidirectional Droplet Propulsion onto Gradient Brushes Without External
Energy Supply

Russell Kajouri  [  Panagiotis E. Theodorakis  [<EMAIL_ADDRESS>Piotr
Deuar  [  Rachid Bennacer  [  Jan Židek  Central European Institute of
Technology, Brno University of Technology, Purkyňova 656/123, 612 00 Brno,
Czech Republic  Sergei A. Egorov  [  Andrey Milchev  Bulgarian Academy of
Sciences, Institute of Physical Chemistry, 1113 Sofia, Bulgaria

###### Abstract

Using extensive molecular dynamics simulation of a coarse-grained model, we
demonstrate the possibility of sustained unidirectional motion (durotaxis) of
droplets without external energy supply when placed on a polymer brush
substrate with stiffness gradient in a certain direction. The governing key
parameters for the specific substrate design studied, which determine the
durotaxis efficiency, are found to be the grafting density of the brush and
the droplet adhesion to the brush surface, whereas the strength of the
stiffness gradient, the viscosity of the droplet or the length of the polymer
chains of the brush have only minor effect on the process. It is shown that
this durotaxial motion is driven by the steady increase of the interfacial
energy between droplet and brush as the droplet moves from softer to stiffer
parts of the substrate whereby the mean driving force gradually declines with
decreasing roughness of the brush surface. We anticipate that our findings
indicate further possibilities in the area of nanoscale motion without
external energy supply.

###### keywords:

Droplets, Gradient Substrates, Durotaxis, Polymer Brush, Motion steering,
Molecular Dynamics

ifpan] Institute of Physics, Polish Academy of Sciences, Al. Lotników 32/46,
02-668 Warsaw, Poland ifpan] Institute of Physics, Polish Academy of Sciences,
Al. Lotników 32/46, 02-668 Warsaw, Poland ifpan] Institute of Physics, Polish
Academy of Sciences, Al. Lotników 32/46, 02-668 Warsaw, Poland ENS]Université
Paris-Saclay, ENS Paris-Saclay, CNRS, LMPS, 4 Av. des Sciences, 91190 Gif-sur-
Yvette, France University Virginia] Department of Chemistry, University of
Virginia, Charlottesville, VA 22901, USA [JoGu University] Institut für
Physik, Johannes Gutenberg Universität Mainz, 55099 Mainz, Germany [Leibniz-
Institut fuer polymerforschung] Leibniz-Institut für Polymerforschung,
Institut Theorie der Polymere, Hohe Str. 6, 01069 Dresden, Germany

##  1 INTRODUCTION

The motion of nano-objects, for example, liquid nanodroplets, can be provoked
and sustained on solid substrates without an external energy supply. Moreover,
the direction of motion can be controlled and nanodroplets can move along
predetermined trajectories. A way of achieving this effect is by placing the
droplet onto a gradient substrate, that is, a substrate with a steadily
varying property along a specific direction. This is particularly attractive
for the development of various technologies in microfluidics,
microfabrication, coatings, nanoscale actuation and energy conversion, and
biology 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12. Various possibilities for the
design of gradient substrates have been reported. For example, durotaxis
motion is caused by changes in stiffness along a substrate, as has been shown
in various natural processes in biology (e.g., cell movement on tissues) 11,
12 and in the case of real and in silico experiments with liquid droplets 13,
14, 15, 16, 17, 18, 19, 20, 21, 22, 23. Another characteristic example is the
rugotaxis motion of droplets on wavy substrates with a gradient in the
wavelength that characterises their pattern 24, 25. Other possibilities
include the use of wettability differences 26, 27 and physical pinning 28.
Recent work has also highlighted the possibility of uni-directional transport
of small condensate droplets on asymmetric pillars29 or three-dimensional
capillary ratchets30. In the latter case, the motion can take place in one or
the other direction, depending on the surface tension of the liquid. Other
possibilities of directional motion can take advantage of charge gradients
that can achieve long-range transport and are based on electrostatic31, 32 or
triboelectric charges.33 In contrast, motion caused by temperature gradient
(thermotaxis) 34, electrical current 35, 36, 37, 38, charge 39, 40, 41, or
even simple stretch 42, would require external energy supply 19, as, also, in
the case of chemically driven droplets 43, 44, droplets on vibrated substrates
45, 46, 47, 48 or wettability ratchets 49, 50, 51, 52.

Inspired by our previous work with specific substrate designs that lead to the
durotaxis13 and rugotaxis24 motion of nanodroplets as motivated by the
corresponding experiments 15, 25, here, we propose a new design for the
substrate, which is capable of sustaining the droplet motion. We consider a
polymer brush, consisting of polymer chains grafted onto a flat, solid
surface, and the stiffness gradient is introduced to the brush substrate by
varying the stiffness of the polymer chains, which in practice amounts to
tuning their persistence length. To understand the mechanism of the durotaxis
motion on brush substrates and analyse the influence of relevant parameters
for the brush and the droplet (e.g., droplet adhesion to the substrate,
droplet size, viscosity, etc.), we have carried out extensive molecular
dynamics (MD) simulations of a coarse-grained (CG) model. This is crucial as
the nanoscale motion of nano-objects is usually controlled by tiny effects at
the interface between the droplet and the substrate resulting from the
molecular interactions between the two, which only a method with molecular
scale resolution can capture. As in the case of durotaxis13 and rugotaxis24
motions, we find that the motion is caused by a gradient in the
droplet–substrate interfacial energy, which translates into an effective force
that drives the droplet towards the stiffer, flatter parts of the substrate.
Moreover, we find that the efficiency of the durotaxis motion for brush
substrates is higher for moderate values of the grafting density and droplet
adhesion to the substrate as well as for smaller droplets and longer brush
chains. Surprisingly, we have not observed a significant effect of the
stiffness gradient in the case of brush substrates when the motion was
successful. We anticipate that our study will shed some light into the
durotaxis motion of droplets on brush, gradient substrates, thus providing
further possibilities in nanoscale science and technology 20, for various
medicine and engineering applications 53, 12. Moreover, brush substrates share
connection with various biological surfaces that expel various exogenous
substances from their structure,54 such as the the mucus layer from airway
epithelia,55 while the gradient concept plays an important role in
applications of regenerative medicine.12 In the following, we discuss our
simulation model and methodology. Then, we will present and discuss our
results, and in the final section we will draw our conclusions.

##  2 MATERIALS AND METHODS

Our system consists of a polymer-brush substrate and a droplet placed on its
soft part (Fig. 1). We have found that the durotaxis motion in the case of
brush substrates takes place from the softer towards its stiffer parts, which
is in line with previous simulation findings for another substrate design with
stiffness gradient.13 In the direction of the stiffness gradient, the
substrate has length, $L_{x}=100~{}\sigma$ ($\sigma$ is the unit of length),
while in the $y$-direction, $L_{y}=50~{}\sigma$, which guarantees that mirror
images of the droplet in this direction will not interact during the course of
the simulation due to the presence of periodic boundary conditions that are
applied in all Cartesian directions. Finally, two walls are placed normal to
the $x$ direction as shown in Fig. 1 and the size of the box in the $x$
direction is large enough to guarantee that there are no interactions between
the walls or the polymers on the two opposite sides of the simulation domain
in the $x$ direction. Wall beads were kept immobile during the simulation.

Figure 1:  Typical initial configuration of the system, where the droplet is
placed on the softest end of the brush substrate. Here, $N_{\rm b}=30$,
$N=4000$, and $N_{\rm d}=10$ beads, $\sigma_{\rm g}=0.6~{}\sigma^{-2}$, and
$\varepsilon_{\rm db}=0.6~{}\epsilon$. At the softest end
$k_{\rm\theta}=0~{}\epsilon/{\rm rad}^{2}$, while at the stiffest end
$k_{\rm\theta}=80~{}\epsilon/{\rm rad}^{2}$, with a linear gradient of
$k_{\theta}$ between them in the $x$ direction. $L_{x}$, $L_{y}$, and $L_{z}$,
indicate the dimensions of the immobile walls. See text for further details.
The snapshot of the system was obtained using Ovito software.56

The standard bead–spring model 57 has been employed in our simulations. In
this model, interactions between any of the system components, i.e. the drop
(d), the brush (b), and the wall (w) beads, are described by means of the
Lennard-Jones (LJ) potential

$U_{\rm LJ}(r)=4\varepsilon_{\rm ij}\left[\left(\frac{\sigma_{\rm
ij}}{r}\right)^{12}-\left(\frac{\sigma_{\rm ij}}{r}\right)^{6}\right],$ (1)

where $r$ is the distance between any pair of beads in the system within a
cutoff distance. Indices ${\rm i}$ and ${\rm j}$ in Eq. 1 indicate the type of
beads. The size of the beads is $\sigma_{\rm ij}=\sigma$ for all interactions.
The LJ potential is cut and shifted at the cutoff distance, $r_{\rm
c}=2.5~{}\sigma$, for the interaction between the droplet (d) beads, as well
as the interaction between the droplet and the brush (b) beads. In contrast, a
purely repulsive potential for the interaction between the brush beads, as
well as between the brush and the wall (w) beads, was considered, that is, in
this case, $r_{\rm c}=2^{1/6}~{}\sigma$. The strength of the attractive
interactions is determined by the parameter $\varepsilon_{\rm ij}$ of the LJ
potential 58 . In our study, $\varepsilon_{\rm dd}=1.5~{}\epsilon$, with
$\epsilon$ defining the energy scale. Moreover, $\varepsilon_{\rm
bb}=\epsilon$, while $\varepsilon_{\rm db}$ is the parameter that controls the
attraction (adhesion) of the droplet to the substrate beads and in our study
ranged from $0.1$ to $1.2~{}\epsilon$.

The grafting density, $\sigma_{\rm g}$, is varied from $0.1$ to
$1.0~{}\sigma^{-2}$ in our study. The size of the droplets can also vary
through the total number of beads that the droplet contains, which ranged
between $2\times 10^{3}$ and $1.6\times 10^{4}$ beads in our simulations.
These beads belong to fully flexible, linear polymer chains. By varying the
length of the droplet chains, e.g., from $10$ to $80$ beads, we can alter the
viscosity of the droplet.13 The finite extensible nonlinear elastic (FENE)
potential 57 was used to tether together consecutive beads in these polymer
chains, as well as the polymer beads along the linear brush-polymer chains.
The mathematical expression for the FENE potential is as follows:

$U_{\rm FENE}(r)=-0.5K_{\rm FENE}R_{\rm 0}^{2}\ln\left[1-\left(\frac{r}{R_{\rm
0}}\right)^{2}\right],$ (2)

where $r$ is the distance between two consecutive beads along the polymer
backbone, while $R_{\rm 0}=1.5~{}\sigma$ expresses the maximum extension of
the bond, and $K_{\rm FENE}=30~{}\epsilon/\sigma^{2}$ is an elastic constant.
Lengths of the polymer chains in the droplet greater than $N_{\rm d}=10$
guarantee that there are no evaporation effects and the vapour pressure is
hence sufficiently low 59 . We have also investigated the effect of the
length, $N_{\rm b}$, of the polymer chains of the brush on the durotaxis
motion, by choosing different lengths, namely $N_{\rm b}=15$, $30$, and $50$
beads.

The stiffness gradient is imposed on the brush substrate by varying the
stiffness of the individual brush polymer-chains. The total length of the
brush chains, $N_{\rm b}$, was the same for all chains, but their stiffness
changed depending on the Cartesian coordinate of their grafting site in the
$x$ direction, i.e. chains with the same position $X$ of their grafted end
have the same stiffness. The chain stiffness was controlled by using a
harmonic angle potential for every triad of consecutive beads along the
polymer chain and tuning its strength through the harmonic constant
$k_{\theta}$. The form of the harmonic potential reads:

$U_{\theta_{\rm ijk}}(\theta)=k_{\theta}(\theta_{\rm ijk}-\theta_{0})^{2},$
(3)

where $\theta_{\rm ijk}$ is the angle between three consecutive beads $\rm i$,
$\rm j$, and $\rm k$ along a brush polymer chain and, $\theta_{0}=\pi~{}{\rm
rad}$, is the equilibrium angle. A linear gradient in the stiffness constant,
$k_{\theta}$, is considered in our study to explore the properties of our
systems. As we will discuss later, while the gradient in the stiffness of the
substrate is necessary to initiate and maintain the durotaxis motion, the
system is rather insensitive to the exact value of the gradient and the key
parameters for the brush substrate turn out to be the grafting density,
$\sigma_{\rm g}$, and the substrate wettability as controlled via the
parameter $\varepsilon_{\rm db}$. The reasons for this will be revealed during
the discussion of our results.

To evolve our system in time, the Langevin thermostat was used, whose details
have been discussed in previous studies60, 61 . Hence, the simulations are in
practice realised in the canonical ensemble,62 where the temperature, $T$, of
the system fluctuates around a predefined value $T=\epsilon/k_{B}$, with
$k_{B}$ being the Boltzmann constant and $\epsilon$ the energy unit. For the
integration of the equations of motion, the LAMMPS package 63 was employed.
The MD time unit is $\tau=\sqrt{m\sigma^{2}/\epsilon}$, where $m$ is the unit
of mass, and the integration time step was $\Delta t=0.005~{}\tau$. Typical
simulation trajectories start from configurations like the one presented in
Fig. 1 with the total length of each trajectory being $10^{8}$ MD integration
steps. If a droplet fully transverses the substrate from the softest to the
stiffest end of the substrate, then the durotaxis motion is considered as
successful. To ensure reliable statistics an ensemble of ten independent
trajectories with different initial conditions (by changing the initial
velocities assigned to each particle) was used for each set of system
parameters. Our results are based on the analysis of these trajectories for
each case.

##  3 RESULTS AND DISCUSSION

Figure 2:  (a) Regime map indicating the probability, $P$ (color scale), that
a droplet will cover the full distance over the substrate in the $x$ direction
from the softest to the stiffest part (successful durotaxis cases) for
different values of the droplet–substrate attraction, $\varepsilon_{\rm db}$,
and the grafting density, $\sigma_{\rm g}$. Probabilities, $P$, are based on
an ensemble of ten independent simulations for each set of parameters. The
regimes where the droplet penetrates into the brush or detaches from the
substrate due to the weak $\varepsilon_{\rm db}$ attraction are also shown
with a different color. (b) The color map indicates the average velocity of
the droplet, $\bar{\varv}=L^{\prime}_{x}/t$, for the successful durotaxis
cases, where $t$ is the time that the droplet needs to cross the full length
of the brush substrate in the $x$ direction, and $L^{\prime}_{x}$ is the
actual distance covered by the centre-of-mass of the droplet for each
successful case. $N=4000$, $N_{\rm d}=10$, $N_{\rm b}=30$ beads. The stiffness
constant for the polymer chains in the softest part of the substrate is zero
(fully flexible chains), growing linearly to $k_{\rm\theta}=80~{}\epsilon/{\rm
rad}^{2}$ at the stiffest part of the substrate. Since $L_{x}=100~{}\sigma$,
the stiffness gradient is $\Gamma=0.8~{}\epsilon/{\rm rad}^{2}\sigma$.
Snapshots on top of the plot indicate examples of detachment, durotaxis, and
penetration, as indicated.  Figure 3:  Density profile in the $z$ direction
for three polymer brushes at $\sigma_{g}=0.6~{}\sigma$ with fully flexible
chains of length, $N_{\rm b}=30$, $60$, and $90$ beads, as indicated. The
inflection point of the curves, $d^{2}\rho/dz^{2}=0$, marked with
$\boldsymbol{\times}$, is shown for each case, which corresponds to the height
of the brush. The position of the centre of mass of the droplet in the $z$
direction is marked with arrows of the same colour for each case. Here,
$N=4000$ and $N_{\rm d}=10$ beads, and $\varepsilon_{\rm db}=0.6~{}\epsilon$.
Snapshots for each case are shown in the plot. In the case of brushes with
longer chains, the droplet is immersed deeper into the brush.

By exploring a wide range of parameters, we have found that the grafting
density, $\sigma_{\rm g}$, and the attraction strength between the droplet and
the substrate, $\varepsilon_{\rm db}$, are two key parameters of the substrate
design, since they greatly affect the possibility for successful durotaxis.
Figure 2 presents the regime maps as a function of these two parameters with
the probability, $P$, of successful durotaxis calculated from the ensemble of
ten independent simulations for each set of parameters. However, in our in
silico experiments, apart from durotaxis motion, we have also documented
situations in which the droplet penetrates into the substrate or detaches from
it. Our results indicate that small attraction strengths will lead to droplet
detachment. This is more probable at smaller grafting densities due to fewer
brush–droplet interactions. In contrast large values of the attraction
strength, $\varepsilon_{\rm db}$, can lead to the penetration of the brush
substrate by the droplet, again for smaller values of the grafting density,
$\sigma_{\rm g}$. Moreover, our results suggest that the difficulty of the
droplet to penetrate into the substrate increases roughly linearly with the
grafting density up to $\varepsilon_{\rm db}=\epsilon$, as evidenced by the
linear boundary in the regime maps. In fact, penetration becomes impossible
when $\sigma_{\rm g}\geq 0.9~{}\sigma^{-2}$, since there is not enough space
among the brush beads to accommodate additional droplet beads. Moreover, the
degree of brush penetration also depends on the length, $N_{\rm b}$, of the
brush chains, as shown in Fig. 3. In general, our data suggest that droplets
are immersed deeper in brushes with longer chains, as shown here in the case
of fully flexible polymers and judging by the centre-of-mass of the droplet
with respect to the position of the brush surface, as defined by the
inflection point at the density profile (Fig. 3). A more detailed study on
this effect could potentially reveal more details for droplets immersed in
brush substrates, but this clearly goes beyond the scope of our current study.

When penetration and detachment of the droplet are avoided, then persistent
durotaxis motion is observed with a certain probability, $P$, which depends on
the choice of $\sigma_{\rm g}$, and $\varepsilon_{\rm db}$ (Fig. 2a). From the
results of Fig. 2a, we find that the range $0.5~{}\sigma^{-2}\leq\sigma_{\rm
g}\leq 0.7~{}\sigma^{-2}$ combined with $0.5~{}\epsilon\leq\varepsilon_{\rm
db}\leq 0.9~{}\epsilon$, in general, provides certainty in the success of the
durotaxis motion ($P=100~{}\%$), since all our droplets were able to fully
cross the substrate from the softest to the stiffest parts of the brush within
the available simulation time of $10^{8}$ MD time steps. As the grafting
density increases above $0.7~{}\sigma^{-2}$, however, we observe that the
probability of durotaxis success suddenly decreases. In this case, a higher
density of the brush chains increases the resultant brush stiffness owing to
the close packing of the chains, which leads to a situation that the role of
the nominal stiffness of the individual brush chains, $k_{\theta}$, in
determining the effective stiffness gradient becomes negligible. As we will
see later in our discussion concerning the underlying durotaxis mechanisms,
the extent of disarray of the brush chain end-monomers at the brush surface is
also reduced pointing to a rather flat density profile.

To determine the efficiency of the durotaxis motion, we have computed the
average velocity, $\bar{\varv}$, of the droplet for the successful cases for
each set of parameters $\sigma_{\rm g}$ and $\varepsilon_{\rm db}$ (Fig. 2b).
Our results indicate that the probability, $P$, for success rather correlates
with the highest average velocity, $\bar{\varv}$, but large values of
$\bar{\varv}$ can also be obtained in certain cases where $P<100~{}\%$, for
example, the case $\sigma_{\rm g}=0.9~{}\sigma^{-2}$, $\varepsilon_{\rm
db}\sim 0.6~{}\epsilon$ (Fig. 2b). This is a clear indication that durotaxis
motion is controlled by tiny effects that can greatly influence the outcome of
the experiments. Moreover, obtaining reliable statistics in cases of $P<100\%$
remains a challenge in MD since this would require the realisation of a large
number of simulations. Hence, as $P\longrightarrow 0$ obtaining reliable
statistics becomes more of a challenge and outliers in the statistics are more
probable. In summary, the plots of Fig. 2 suggest that if one would like all
droplets to fully cross the substrate in the direction of the stiffness
gradient in the shortest time, then values of $\sigma_{\rm
g}=0.6~{}\sigma^{-2}$ and $\varepsilon_{\rm db}\approx 0.7~{}\epsilon$ would
constitute an optimal choice in the in silico experiments. Hence, we argue
that moderate values of $\sigma_{\rm g}$ and $\varepsilon_{\rm db}$ favour
successful and efficient (in terms of time to cross the whole substrate)
durotaxis motion.

Figure 4:  (a) Regime map indicating the probability, $P$ (color scale), that
a droplet will cover the full distance over the substrate in the $x$ direction
from the softest to the stiffest part (sucessful durotaxis cases) for
different values of the droplet–substrate attraction, $\varepsilon_{\rm db}$,
and the stiffness gradient $\Gamma=dk_{\theta}/dx$. The regimes where the
droplet penetrates into the brush or detaches from the substrate due to the
weak $\varepsilon_{\rm db}$ attraction are also shown with a different color.
(b) The color map indicates the average velocity of the droplet,
$\bar{\varv}=L^{\prime}_{x}/t$, for the successful durotaxis cases, where $t$
is the time that the droplet needs to cross the full length of the brush
substrate in the $x$ direction, and $L^{\prime}_{x}$ is the actual distance
covered by the centre-of-mass of the droplet for each successful case. The
stiffness constant, $k_{\theta}$, for the polymer chains in the softest part
of the substrate is zero (fully flexible chains), growing linearly to its
highest value at $x=L_{x}=100~{}\sigma$, which depends on the chosen stiffness
gradient, $\Gamma$. Here, $\sigma_{\rm g}=0.6~{}{\sigma^{-2}}$, $N=4000$,
$N_{\rm b}=30$, and $N_{\rm d}=10$ beads.

Finally, we have explored the effect of the stiffness gradient on the
durotaxis motion. Here, we have picked the best case of Fig. 2, that is
$\sigma_{\rm g}=0.6~{}\sigma^{-2}$ and $\varepsilon_{\rm db}=0.6~{}\epsilon$,
and varied the stiffness gradient $\Gamma=dk_{\theta}/dx$ in the range 0.2–2.6
for fully flexible chains at the softest part of the substrate, which has
actually provided the best result in terms of durotaxis success and
efficiency. Overall, we have found that durotaxis is insensitive to the value
of $\Gamma$ in the range $0.5~{}\epsilon\leq\varepsilon_{\rm db}\leq
0.9~{}\epsilon$ (Fig. 4a), in contrast to what has been observed for other in
silico substrate designs13. Moreover, the average velocities are spread out
with small variations and no indication of a clear trend (Fig. 4b) that would
indicate that a larger stiffness gradient would lead to more efficient
durotaxis motion exists, which has been the case for other substrate
designs13. Moreover, since the motion is most efficient when the softest part
consists of fully flexible chains ($k_{\theta}=0~{}\epsilon/{\rm rad}^{2}$)
might suggest that brush substrates with polymer chains of small persistence
lengths (as soft as possible) are more suitable for successful durotaxis
motion. Henceforth, all our results refer to brush substrates with fully
flexible chains at their softest part and stiffness gradient
$\Gamma=0.8~{}\epsilon/{\rm rad}^{2}\sigma$.

Figure 5:  Average droplet velocity, $\bar{\varv}$, as a function of
$\varepsilon_{\rm db}$ for different (a) droplet size, $N$ ($N_{\rm d}=10$,
$N_{\rm b}=30$ beads), (b) chain length of the polymer chains of the droplet,
$N_{\rm d}$ ($N=4000$ and $N_{\rm b}=30$ beads) , and (c) brush polymer chain
length, $N_{\rm b}$ ($N=4000$ and $N_{\rm d}=10$ beads), as indicated. (d)
Documented maximum average velocity, $\bar{\varv}$, indicated by the colour
map, as a function of $\varepsilon_{\rm db}$ and $N$. Inset shows the average
velocity, $\bar{\varv}$, for all pairs of ($\varepsilon_{\rm db}$, $N$). The
black points indicate the pairs ($\varepsilon_{\rm db}$, $N$) for which we
have the maximum average velocity, $\bar{\varv}$, which is shown in the main
plot. (e) Same as (d), but data are plotted as a function of
($\varepsilon_{\rm db}$, $N_{\rm d}$). (f) In this case data are plotted as a
function of ($\varepsilon_{\rm db}$, $N_{\rm b}$). $\sigma_{\rm
g}=0.6~{}\sigma^{-2}$, $\Gamma=0.8~{}\epsilon/{\rm rad}^{2}\sigma$ in all
cases. Panels (a-c) here show data as a function of $\varepsilon_{\rm db}$
with one other quantity varying, the others held constant. Therefore, it does
take care to carry out a systematic study with as much held constant between
data as useful. Panels (d-f) visualize the same data as in the top panels but
with a different visualization in terms of maximum velocities.

In the following, we examine the effect of various parameters on the
efficiency of the durotaxis motion. For our analysis, we have picked the case
$\sigma_{\rm g}=0.6~{}\sigma^{-2}$, which has shown the best performance in
terms of the probability, $P$, and the average velocity, $\bar{\varv}$, in our
study (Fig. 2) and therefore would most probably allow for exploring a wider
range of the parameter space. Then, Figs 5a,d illustrate the dependence of the
average velocity, $\bar{\varv}$, on the attraction strength, $\varepsilon_{\rm
db}$, for various droplet sizes, $N$. We observe that the fastest durotaxis
motion takes place for the smallest droplets, namely $N=2000$ beads. As
$\varepsilon_{\rm db}$ increases, a maximum value of the average velocity,
$\bar{\varv}$, appears for each case with droplets of different size. In
addition, we consistently see that this maximum value of $\bar{\varv}$ becomes
smaller for larger droplets. For example, for the droplet of $N=2000$ beads
the average velocity measured over ten trajectories is $\bar{\varv}\approx
4.7\times 10^{-4}~{}\sigma/\tau$, while $\bar{\varv}\approx 3.5\times
10^{-4}~{}\sigma/\tau$ for droplets with $16000$ beads. These differences are
generally considered small, especially compared to other in silico
experiments.13, 24 Moreover, we see that the maximum shifts to higher values
of the attraction strength, $\varepsilon_{\rm db}$, as the size of the droplet
increases. For example, the maximum velocity is observed when
$\varepsilon_{\rm db}=0.4~{}\epsilon$ for droplets with $N=2000$ beads and for
$\varepsilon_{\rm db}=0.9~{}\epsilon$ for droplets of $16000$ beads (Fig. 5d).
Finally, for small droplets we observe a steep increase of the average
velocity with $\varepsilon_{\rm db}$, and then a slow decrease (Fig. 3a). For
medium-size droplets, (i.e. $N=4000$ and $N=8000$ beads), there is a smooth
maximum that develops in the middle range of $\varepsilon_{\rm db}$, i.e.
$0.6~{}\epsilon\leq\varepsilon_{\rm db}\leq 0.7~{}\epsilon$, while in the case
of droplets with $N=16000$ beads there is a maximum that slowly develops as
$\varepsilon_{\rm db}$ increases, which is followed by a steeper decrease when
$\varepsilon_{\rm db}>0.9~{}\epsilon$. In summary, we observe that the size of
the droplets is an important parameter for the durotaxis motion.

The next parameter to examine is the chain length, of the polymers comprising
the droplet, $N_{\rm d}$ (Figs 5b, e). In practice, longer chain lengths would
result in a larger droplet viscosity. During this and previous work with
polymer liquid-droplets13 we have determined that the most relevant values for
our study are within the range $10\leq N_{\rm d}\leq 80$ beads. Interestingly,
the droplet viscosity seems not to play an important role in the overall
efficiency of the durotaxis motion, for a given attraction strength,
$\varepsilon_{\rm db}$. In other words, droplets with different $N_{\rm d}$
would exhibit a similar durotaxial efficiency for a specific choice of
$\varepsilon_{\rm db}$. As a result, the maximum average velocity,
$\bar{\varv}$, as a function of $\varepsilon_{\rm db}$ appears at
$\varepsilon_{\rm db}\sim 0.6-0.7~{}\epsilon$ (Fig. 5e). Hence, we can
conclude that droplets with different viscosity will have a similar durotaxis
performance and, here, a moderate choice for the value of droplet–substrate
attraction strength would yield the fastest durotaxis motion.

The effect of the length of the brush chains, $N_{\rm b}$, on the durotaxis
motion is shown in Figs 5c, f. We find that the larger the length, $N_{\rm
b}$, the more efficient the durotaxis motion becomes. Although a larger
difference is noticed when $N_{\rm b}$ was doubled from 15 to 30 beads, a
saturation in our data occurred when $N_{\rm b}$ increased from 30 to 50
beads. Overall, our results indicate that brushes with longer polymer chains
favor the durotaxis motion. Moreover, for $N_{\rm b}=30$ and 50 beads, the
maximum velocity is found for $\varepsilon_{\rm db}=0.6~{}\epsilon$, and hence
is independent of the choice of $N_{\rm b}$. The reasons for this behavior
will become more apparent as we will discuss the underlying mechanism of the
droplet motion in the following.

Figure 6:  Droplet–brush interfacial energy, $U_{\rm db}$, as a function of
time, $t$ (a, dashed line is a guide for the eye), and the position, $X$, of
the centre of mass of the droplet in the $x$ direction (b), for a case with
successful durotaxis ($N=4000$, $N_{\rm d}=10$, and $N_{\rm d}=30$ beads.
$\varepsilon_{\rm db}=0.6~{}\epsilon$ and $\sigma_{\rm g}=0.6~{}\sigma^{-2}$).
Inset shows the force, $f_{x}=-\frac{\partial U_{\rm db}}{\partial x}$ based
on a nonlinear fit of the $\tanh$ function on the $U_{\rm db}$ data. The
higher concentration of points in panel (b) at certain ranges of $X$ simply
indicates that the droplet spends more time at these positions as it moves to
the stiffest parts of the substrate. The fit function of the $U_{\rm db}$ only
provides an average picture of the decay of the interfacial energy.  Figure 7:
(a) The standard deviation in the end-to-end distance of the brush polymer
chains as a function of their grafting position, $X$. Larger values of $X$
correspond to the stiffer parts of the substrate. Inset shows the density
profile in the $z$ direction at different positions, $X$. The colour reflects
the stiffness of the chains ($k_{\theta}$). $\varepsilon_{\rm
db}=0.6~{}\epsilon$, $\sigma_{\rm g}=0.6~{}\sigma^{-2}$.
$\Gamma=0.8~{}\epsilon/{\rm rad}^{2}\sigma$ and initial stiffness
$0~{}\epsilon/{\rm rad}^{2}$ at the softest end. $N=4000$, $N_{\rm d}=10$,
$N_{\rm b}=30$ beads. (b) Interpenetration length, $W$. $l_{z}$ is the average
distance of contact pairs between the droplet and the brush for the same set
of parameters as in (a) with $X$ here indicating the centre-of-mass position
of the droplet in the $x$ direction along the gradient. (c) Average
interpenetration length as a function of the grafting density, $\sigma_{\rm
g}$, for substrates with constant stiffness ($k_{\theta}=0~{}\epsilon/{\rm
rad}^{2}$ or $k_{\theta}=80~{}\epsilon/{\rm rad}^{2}$) or with stiffness
gradient $\Gamma=0.8~{}\epsilon/{\rm rad}^{2}\sigma$ and initial stiffness
$0~{}\epsilon/{\rm rad}^{2}$ at the softest end, as indicated. $N=4000$,
$N_{\rm d}=10$, $N_{\rm b}=30$ beads. $\varepsilon_{\rm db}=0.6~{}\epsilon$ in
the main plot and $\sigma_{\rm g}=0.6~{}\sigma^{-2}$ in the inset.

From our earlier studies, we have seen that the minimization of the
interfacial energy, $U_{\rm db}$, between the droplet and the substrate is the
driving force for the durotaxis motion in the case of substrates with
stiffness gradient 13 or wrinkle substrates with a gradient in the wavelength
characterising the wrinkles 24. This has also been found in the case of a
nanoflake on substrates with stiffness gradient.16 Moreover, we have argued
that the efficiency of the motion depends on the rate of change of the
interfacial energy along the gradient direction. Hence, it is relevant to
examine the interfacial energy, $U_{\rm db}$, as a function of time, $t$, and
the coordinate, $X$, of the centre-of-mass of the droplet along the substrate
(see Fig. 6 showing these quantities for a successful durotaxis case). Our
results indicate that $U_{\rm db}$ over time reaches more negative values
(Fig. 6a), which corresponds to a larger number of attractive pair
interactions between the droplet and the substrate. Importantly, we also see
that the energy decreases as a function of the position of the droplet, $X$,
which clearly indicates that the droplet moves to areas of more negative
energy (Fig. 6b) towards the stiffer parts of the substrate (increasing $X$).
A larger decrease of the energy takes place in the initial soft parts of the
substrate at small $X$ values, while the change in the substrate–droplet
interfacial energy is smaller after the droplet has moved a distance about 40
$\sigma$. Moreover, we can see that the droplet does not spend the same time
at each position $X$, as seen by the density of points at certain positions,
$X$, and knowing that samples have been taken over time at equal intervals.
Hence, the durotaxis motion of the droplet cannot be characterised as a steady
state. To illustrate the effect of the gradient in the interfacial energy,
$U_{\rm db}$, one can actually plot the negative derivative of the interfacial
energy (inset of Fig. 6b) after performing a suitable fit on the $U_{\rm db}$.
This derivative would correspond to an on-average driving force
$f_{x}=-\frac{\partial U_{\rm db}}{\partial x}$ that propels the droplet
toward the stiffer parts of the brush substrate. It is clear that this ‘force’
becomes significantly smaller when the droplet enters the stiffer half of the
substrate in the $x$ direction. However, this small ‘force’ is still able to
sustain the motion of the droplet and lead to a successful durotaxis case.
Still, the inset of Fig. 6b reflects the observed behavior caused by an
average ‘force’ that propels the droplet. Finally, we have clearly seen from
our data that unsuccessful durotaxis cases are characterised by a flat
interfacial energy that fluctuates around a constant value, which would
ideally yield $f_{\rm x}\approx 0$.

To further understand the durotaxis mechanism of droplets on brush substrates,
we have gone one step further and tried to identify the origin of the changing
interfacial energy in successful durotaxis cases. In particular, we have
measured the standard deviation of the end-to-end distance, $\Omega$, which
describes the width of the free-end positions distribution of the brush
polymer chains (Fig. 7a), i.e., the surface roughness of the brush. We can
observe that $\Omega$ decays monotonically towards the stiffest parts of the
brush. In other words, stiffer polymer chains exhibit a smaller extent of
fluctuations concerning the end-to-end distance of the polymer chains, which
is generally expected when the stiffness increases. Moreover, we can see that
the decay of $\Omega$ is faster at the soft parts of the substrate and
generally follows the decay in the interfacial energy (Fig. 6). The inset of
Fig. 7a presents results for the the density profile in the $z$ direction at
different positions $X$ in the $x$ direction (along the stiffness gradient).
We observe that the thickness of the brush surface becomes smaller towards the
stiffer parts, which would correspond to a flatter surface locally. In
contrast, slower decaying density profiles correspond to a larger thickness of
the brush surface, in other words, to a rougher brush surface, which is
observed in the soft substrate parts. In practice, rough surfaces result in a
smaller number of contacts with the droplet, and as a result a higher (less
negative) interfacial energy.13 On the contrary, a flat profile would allow
for a larger number of contacts between the droplet and the substrate. For
this reason, the droplet moves towards the stiffer parts of the substrate. We
provide further evidence for our argument by measuring the interpenetration
length, $W$, as a function of the position, $X$, of the centre-of-mass of the
droplet along the substrate in the direction of the stiffness gradient (Fig.
7b). This property reflects the average distance of the substrate–droplet
contact pairs, which is noted here with the symbol $l_{z}$. We can clearly see
that the interpenetration length decreases at the stiffer parts of the
substrate, which points to a sharper (flatter) brush surface, in accordance
with the results of Fig. 7a. We have further explored the dependence of $W$ on
$\varepsilon_{\rm db}$ and $\sigma_{g}$ and found that it decreases as a
function of $\sigma_{\rm g}$ due to the induced stiffness by the steric
interactions between the brush polymer chains, while it increases with
$\varepsilon_{\rm db}$. An important conclusion from the results of Fig. 7c is
that $W$ is clearly larger in the case of soft brushes
($k_{\theta}=0~{}\epsilon/{\rm rad}^{2}$). However, it is the gradual change
of this roughness that plays an important role in inducing the interfacial
gradient, which in turn translates into the effective ‘force’ that drives the
droplet motion.

Figure 8:  Droplet–brush interfacial energy, $U_{\rm db}$, per area of the
droplet–substrate contact surface, $A_{\rm db}$, as a function of
$\varepsilon_{\rm db}$. $\sigma_{\rm g}=0.6~{}\sigma^{-2}$ and
$\Gamma=0.8~{}\epsilon/{\rm rad}^{2}\sigma$. Panels show results for different
(a) $N$ ($N_{\rm b}=30$ and $N_{\rm d}=10$ beads), (b) $N_{\rm d}$ ($N_{\rm
b}=30$ beads), and (c) $N_{\rm b}$ ($N_{\rm d}=10$ beads). The area, $A_{\rm
db}$, was calculated by using the qhull library64, taking the area of the
convex hull.  Figure 9:  (a) Droplet–brush average interfacial energy, $U_{\rm
db}$, per area, $A_{\rm db}$, as a function of the grafting density,
$\sigma_{\rm g}$. Dashed-dotted line is the result of a linear fit with slope
$-0.57~{}\epsilon$ as indicated. Inset illustrates the values of the area,
$A_{\rm db}$, versus the grafting density. $\varepsilon_{\rm
db}=0.6~{}\epsilon$, $N=4000$, $N_{\rm b}=30$ and $N_{\rm d}=10$ beads.
$\Gamma=0.8~{}\epsilon/{\rm rad}^{2}\sigma$. (b) Advancing, $\theta_{\rm A}$,
and receding, $\theta_{\rm R}$, contact angles along the durotaxis motion ($x$
direction) as a function of $\varepsilon_{\rm db}$ for droplets of different
size as indicated. $\sigma_{\rm g}=0.6~{}\sigma^{-2}$. The inset shows the
dependence of the contact angles on the grafting density, $\sigma_{\rm g}$,
for $\varepsilon_{\rm db}=0.6~{}\epsilon$. In both the main panel and the
inset, $N_{\rm d}=10$ and $N_{\rm b}=30$ beads, and
$\Gamma=0.8~{}\epsilon/{\rm rad}^{2}\sigma$.

With these indications that the minimization of the interfacial energy drives
the droplet toward areas of more negative energy, in Fig. 8, we show the
interfacial energy for the same systems shown in Fig. 5. When the chain length
of the droplet polymer chains, $N_{\rm d}$, varies within the range considered
in our study, there are no noticeable changes in the interfacial energy (Fig.
8b) and systems are rather characterised by an equivalent average energy
profile, which is reflected in the average velocity of the droplets (Fig. 5b).
There, we have found that a change in the viscosity of the droplet does not
significantly affect the durotaxis efficiency. Slight changes in the average
interfacial energy are observed when the brush chain length, $N_{\rm b}$,
varies (Fig. 8c). These results are also consistent with those of Fig. 5c. In
particular, we can see that the brush with chain length $N_{\rm b}=50$ beads
appears to have the lowest average interfacial energy and the largest
interfacial area for the whole range of $\varepsilon_{\rm db}$ considered
here. Moreover, differences in the interfacial energy between the systems with
different $N_{\rm b}$ become smaller as $\varepsilon_{\rm db}$ increases,
which is much clearer in the data concerning the interfacial area, $A_{\rm
db}$, in line with the results of Fig. 5c. In contrast, rather larger
differences are observed when the size of the droplets, $N$, changes (Fig.
8a). In particular, droplets of smaller size have a markedly lower (less
negative) interfacial energy per area than larger droplets. These differences
become more apparent as the strength of interactions, $\varepsilon_{\rm db}$,
increases. Interestingly, although the energy here shows a monotonic behavior,
this is not the case for the average velocity of the droplet (Fig. 5a), which
indicates that droplets of different size are differently affected in a global
sense for a given set of substrate parameters. One contributing factor is that
the change in interfacial energy from soft to hard substrate is proportional
to the contact area $\propto N^{2/3}$, while the mass is proportional to $N$.
Therefore the energy change per atom is $\propto N^{-1/3}$, other factors
being equal, and the velocity change can be expected to scale as $N^{-1/6}$.
This is not inconsistent with Fig. 5a at medium strength $\varepsilon_{\rm
db}$. Finally, the interfacial energy decreases proportionally to the grafting
density, with a slope of $-0.57~{}\epsilon$ (Fig. 9a), suggesting a
proportionally larger number of contacts between the droplet and the substrate
as the grafting density grows.

We have also monitored the advancing, $\theta_{\rm A}$, and receding,
$\theta_{\rm R}$, contact angles during the durotaxial motion and results are
presented in Fig. 9b for typical cases. In this case, the angles have been
determined by using the curvature of the droplet as described in previous
studies65, thus avoiding error-prone fits. We observe that both $\theta_{\rm
A}$ and $\theta_{\rm R}$ decrease rather linearly with the increase of
$\varepsilon_{\rm db}$. A linear decrease has been observed for droplets on
solid substrates.65, 13 Droplets of different size, namely $N=2000$ and
$N=16000$ beads show the same trend. Moreover, as $\varepsilon_{\rm db}$
increases, the differences between smaller and larger droplets become less
pronounced. A much weaker dependence of the contact angles on the grafting
density has been observed for $\sigma_{\rm g}\leq 0.6~{}\sigma^{-2}$ and a
small scale rather linear dependence for $\sigma_{\rm g}\geq
0.6~{}\sigma^{-2}$. Our results, which are averaged over the whole trajectory,
do not show any statistically significant difference between advancing and
receding contact angles. Due to the greater magnitude of interfacial energy
(more negative) at the stiffer parts of the substrate (Fig. 6) one would
expect a smaller advancing contact angle (for example, as Fig. 9b suggests,
larger attraction leads to smaller contact angles). Slightly larger values for
the receding contact angle, $\theta_{\rm R}$, are observed consistently for
different values of $\sigma_{\rm g}$ and $\varepsilon_{\rm db}$, but, still,
within the statistical error. In experiment and theory, when the droplet moves
by steadily applied external force the receding contact angles are smaller
than the advancing ones in which case friction effects might also play a role.

Figure 10:  Diffusion coefficient ($\sigma_{\rm g}=0.6~{}\sigma^{-2}$) as a
function of (a) $\varepsilon_{\rm db}$ (the inset shows the dependence on the
grafting density, $\sigma_{\rm g}$ for the case $\varepsilon_{\rm
db}=0.6~{}\epsilon$). $N_{\rm d}=10$ and $N_{\rm b}=30$ beads. (b) Diffusion
coefficient for different droplet chain length, $N_{\rm d}$, brush chain
length, $N_{\rm b}$ (inset) as a function of $\varepsilon_{\rm db}$.
$\sigma_{\rm g}=0.6~{}\sigma^{-2}$, $\varepsilon_{\rm db}=0.6~{}\epsilon$ in
the inset. $N=4000$ beads. (c) Typical, successful durotaxis trajectories of
the centre of mass of the droplet in the $x-y$ space as indicated by different
colours for separate run. $\sigma_{\rm g}=0.6~{}\sigma^{-2}$,
$\varepsilon_{\rm db}=0.6~{}\epsilon$, $N=4000$, $N_{\rm b}=30$ and $N_{\rm
d}=10$ beads, and $\Gamma=0.8~{}\epsilon/{\rm rad}^{2}\sigma$.

Finally, to investigate the diffusion of the droplet, and how it is affected
by the main parameters, we have analysed the movement of the centre of mass of
the droplet in the $y$ direction (Fig. 10). Dynamics in the $y$ direction is
free of extraneous effects due to grafting density gradient and the durotaxis
itself, so serves as a good control and measurement of the diffusion
properties. Overall droplets with higher viscosity and larger wettability on
the brush substrate, show a slower diffusion. The same effect has the grafting
density. In contrast, $N_{\rm b}$ does not have a tangible effect on the
droplet motion. In Fig. 10c, we present typical trajectories as a function of
the $X$ and $Y$ positions of the droplet centre-of-mass. We can observe that
durotaxis motion is also characterised by significant random motion in both
$x$ and $y$ axes including ‘reversals’. In some cases, we can observe that the
droplet can cover in the $y$ direction a distance as much as 30% of the total
distance in the $x$ direction, which suggests that durotaxis motion is also
strongly affected by random effects at the interface between the droplet and
the substrate, especially when the gradient of the interfacial energy becomes
small, i.e, the driving force of the durotaxis is correspondingly small.

##  4 CONCLUSIONS

In this study, we have proposed a new design of brush polymer substrates that
is capable of leading to the durotaxis motion of nanodroplets. The knowledge
gained here may lead to new experimental brush-based substrate designs and
provide further understanding of relevant biological processes, such as the
motion of cells on tissues11, 12, or the mucus flow around lung cilia55. Our
analysis has also indicated that the durotaxis motion on brush substrates is
driven by a corresponding gradient in the interfacial energy between the
droplet and the substrate, in line with previous findings in the context of
various other substrate designs.13, 24, 16 Moreover, we have found that the
origin of the steady increase of the interfacial energy is related to the
state of the brush surface, which appears as a ‘rougher’ profile in the softer
parts of the substrate and a flatter interface in the stiffer regions. This
translates into a larger number of contacts between the droplet and the
substrate in the stiffer parts, and hence a more negative interfacial energy
along the direction of the stiffness gradient.

We have also conducted a parametric study based on the various system
parameters in order to gain further insight into the system and identify the
key parameters of the brush substrate design. Two of the key parameters are
the grafting density of the polymer chains and the substrate wettability. Our
findings suggest that the durotaxis motion is favored by a moderate grafting
density, which in our case translates to values of grafting density
$\sigma_{\rm g}\approx 0.5~{}\sigma^{-2}$, and the same time moderate values
of the substrate wettability, namely $0.5~{}\epsilon\leq\varepsilon_{\rm
db}\leq 0.9~{}\epsilon$. Surprisingly, we have found that the stiffness
gradient itself, as defined by the linear change in the stiffness of
individual polymer chains tuned by the harmonic constant, $k_{\rm\theta}$,
does not induce important changes in the efficiency of the durotaxis motion.
This might be due to the relatively large grafting density, which affects the
apparent stiffness of the polymer chains, thus minimising the actual effect of
the gradient. Such an effect has not been observed in previous systems.13
Moreover, the droplet viscosity also seems not to affect the durotaxis
efficiency as also the length of the polymer brush chains, since the durotaxis
performance seems to saturate after a certain polymer length. In contrast, the
size of the droplet plays a role. In particular, smaller droplets seem to
reach a faster durotaxis and for a lower adhesion to the subtrate, while
larger droplets move with a lower average velocity exhibiting their maximum
velocity at larger droplet–substrate adhesion. Furhtermore, we have not
identified any tangible differences between the advancing and receding contact
angles during the durotaxis motion. The durotaxis motion is induced by tiny
effects at the droplet–substrate interface as judged by the gradient in the
droplet–substrate interfacial energy. Hence, there is a back and forth
wiggling motion while the brush slowly guides the droplet towards a lower
energy state. Finally, we have discussed how our findings could motivate
further experimental research in the area of self-sustained fluid motion on
brush gradient substrates. It would also be interesting to investigate the
droplet durotaxis behavior when many droplets are placed onto the substrate
and explore various effects, such as droplet coalescence or how the
droplet–substrate interactions are affected in populations of droplets.
Further work to explore such possibilities is expected in this direction in
the future. We anticipate that this study for the first time presents new
possibilities for the implementation and understanding of durotaxis motion of
fluids on brush substrates with important implications for various areas, for
example, in the context of biology.

This research has been supported by the National Science Centre, Poland, under
grant No. 2019/35/B/ST3/03426. A. M. acknowledges support by COST (European
Cooperation in Science and Technology [See http://www.cost.eu and
https://www.fni.bg] and its Bulgarian partner FNI/MON under KOST-11). This
research was supported in part by PLGrid Infrastructure.

M*.mp4: Movie illustrates an example of durotaxis motion onto a gradient brush
substrate. The parameters for the system are: $N=4000$, $N_{\rm d}=10$, and
$N_{\rm b}=30$ beads. $\sigma_{\rm g}=0.6~{}\sigma^{-2}$,
$\Gamma=0.8~{}\epsilon/{\rm rad}^{2}\sigma$ (chains at the very soft right end
are fully flexible, i.e., $k_{\theta}=0~{}\epsilon/{\rm rad}^{2}$, while at
the very stiff part $k_{\theta}=80~{}\epsilon/{\rm rad}^{2}$), and
$\varepsilon_{\rm db}=0.6~{}\epsilon$.

## References

  * Srinivasarao et al. 2001 Srinivasarao, M.; Collings, D.; Philips, A.; Patel, S. Three-dimensionally ordered array of air bubbles in a polymer film. _Science_ 2001, _292_ , 79–83
  * Chaudhury and Whitesides 1992 Chaudhury, M. K.; Whitesides, G. M. How to Make Water Run Uphill. _Science_ 1992, _256_ , 1539–1541
  * Wong et al. 2011 Wong, T.-S.; Kang, S. H.; Tang, S. K. Y.; Smythe, E. J.; Hatton, B. D.; Grinthal, A.; Aizenberg, J. Bioinspired self-repairing slippery surfaces with pressure-stable omniphobicity. _Nature_ 2011, _477_ , 443–447
  * Lagubeau et al. 2011 Lagubeau, G.; Le Merrer, M.; Clanet, C.; Quéré, D. Leidenfrost on a ratchet. _Nat. Phys._ 2011, _7_ , 395–398
  * Prakash et al. 2008 Prakash, M.; Quéré, D.; Bush, J. W. Surface tension transport of prey by feeding shorebirds: The capillary ratchet. _Science_ 2008, _320_ , 931–934
  * Darhuber and Troian 2005 Darhuber, A.; Troian, S. Principles of microfluidic actuation by modulation of surface stresses. _Annu. Rev. Fluid Mech._ 2005, _37_ , 425–455
  * Yao and Bowick 2012 Yao, Z.; Bowick, M. J. Self-propulsion of droplets by spatially-varying surface topography. _Soft Matter_ 2012, _8_ , 1142–1145
  * Li et al. 2018 Li, H.; Yan, T.; Fichthorn, K. A.; Yu, S. Dynamic contact angles and mechanisms of motion of water droplets moving on nano-pillared superhydrophobic surfaces: A molecular dynamics simulation study. _Langmuir_ 2018, _34_ , 9917–9926
  * Becton and Wang 2016 Becton, M.; Wang, X. Controlling nanoflake motion using stiffness gradients on hexagonal boron nitride. _RSC Adv._ 2016, _6_ , 51205–51210
  * van den Heuvel and Dekker 2007 van den Heuvel, M. G. L.; Dekker, C. Motor proteins at work for nanotechnology. _Science_ 2007, _317_ , 333–336
  * DuChez et al. 2019 DuChez, B. J.; Doyle, A. D.; Dimitriadis, E. K.; Yamada, K. M. Durotaxis by human cancer cells. _Biophys. J._ 2019, _116_ , 670–683
  * Khang 2015 Khang, G. Evolution of gradient concept for the application of regenerative medicine. _Biosurface Biotribology_ 2015, _1_ , 202–213
  * Theodorakis et al. 2017 Theodorakis, P. E.; Egorov, S. A.; Milchev, A. Stiffness-guided motion of a droplet on a solid substrate. _J. Chem. Phys._ 2017, _146_ , 244705
  * Lo et al. 2000 Lo, C.-M.; Wang, H.-B.; Dembo, M.; Wang, Y.-L. Cell movement is guided by the rigidity of the substrate. _Biophys. J._ 2000, _79_ , 144–152
  * Style et al. 2013 Style, R. W.; Che, Y.; Park, S. J.; Weon, B. M.; Je, J. H.; Hyland, C.; German, G. K.; Power, M. P.; Wilen, L. A.; Wettlaufer, J. S.; Dufresne, E. R. Patterning droplets with durotaxis. _Proc. Natl. Acad. Sci. U.S.A._ 2013, _110_ , 12541–12544
  * Chang et al. 2015 Chang, T.; Zhang, H.; Guo, Z.; Guo, X.; Gao, H. Nanoscale directional motion towards regions of stiffness. _Phys. Rev. Lett._ 2015, _114_ , 015504
  * Pham et al. 2016 Pham, J. T.; Xue, L.; Del Campo, A.; Salierno, M. Guiding cell migration with microscale stiffness patterns and undulated surfaces. _Acta Biomaterialia_ 2016, _38_ , 106–115
  * Lazopoulos and Stamenović 2008 Lazopoulos, K. A.; Stamenović, D. Durotaxis as an elastic stability phenomenon. _J. Biomech._ 2008, _41_ , 1289–1294
  * Becton and Wang 2014 Becton, M.; Wang, X. Thermal gradients on graphene to drive nanoflake motion. _J. Chem. Theory Comput._ 2014, _10_ , 722–730
  * Barnard 2015 Barnard, A. S. Nanoscale locomotion without fuel. _Nature_ 2015, _519_ , 37–38
  * Palaia et al. 2021 Palaia, I.; Paraschiv, A.; Debets, V. E.; Storm, C.; Šarić, A. Durotaxis of Passive Nanoparticles on Elastic Membranes. _ACS Nano_ 2021, _15_ , 15794–15802
  * Tamim and Bostwick 2021 Tamim, S. I.; Bostwick, J. B. Model of spontaneous droplet transport on a soft viscoelastic substrate with nonuniform thickness. _Phys. Rev. E_ 2021, _104_ , 034611
  * Bardall et al. 2020 Bardall, A.; Chen, S.-Y.; Daniels, K. E.; Shearer, M. Gradient-induced droplet motion over soft solids. _IMA J. Appl. Math_ 2020, _85_ , 495–512
  * Theodorakis et al. 2002 Theodorakis, P. E.; Egorov, S. A.; Milchev, A. Rugotaxis: Droplet motion without external energy supply. _EPL_ 2002, _137_ , 43002
  * Hiltl and Böker 2016 Hiltl, S.; Böker, A. Wetting Phenomena on (Gradient) Wrinkle Substrates. _Langmuir_ 2016, _32_ , 8882–8888
  * Pismen and Thiele 2006 Pismen, L. M.; Thiele, U. Asymptotic theory for a moving droplet driven by a wettability gradient. _Phys. Fluids_ 2006, _18_ , 042104
  * Wu et al. 2017 Wu, H.; Zhu, K.; Cao, B.; Zhang, Z.; Wu, B.; Liang, L.; Chai, G.; Liu, A. Smart design of wettability-patterned gradients on substrate-independent coated surfaces to control unidirectional spreading of droplets. _Soft Matter_ 2017, _13_ , 2995–3002
  * Theodorakis et al. 2021 Theodorakis, P. E.; Amirfazli, A.; Hu, B.; Che, Z. Droplet Control Based on Pinning and Substrate Wettability. _Langmuir_ 2021, _37_ , 4248–4255
  * Feng et al. 2020 Feng, S.; Delannoy, J.; Malod, A.; Zheng, H.; Quéré, D.; Wang, Z. Tip-induced flipping of droplets on Janus pillars: From local reconfiguration to global transport. _Sci. Adv._ 2020, _6_ , eabb5440
  * Feng et al. 2021 Feng, S.; Zhu, P.; Zheng, H.; Zhan, H.; Chen, C.; Li, J.; Wang, L.; Yao, X.; Liu, Y.; Wang, Z. Three dimensional capillary ratchet-induced liquid directional steering. _Science_ 2021, _373_ , 1344–1348
  * Sun et al. 2019 Sun, Q.; Wang, D.; Li, Y.; Zhang, J.; Ye, S.; Cui, J.; Chen, L.; Wang, Z.; Butt, H. J.; Vollmer, D.; Deng, X. Surface charge printing for programmed droplet transport. _Nat. Mater._ 2019, _18_ , 936–941
  * Jin et al. 2022 Jin, Y.; Xu, W.; Zhang, H.; Li, R.; Sun, J.; Yang, S.; Liu, M.; Mao, H.; Wang, Z. Electrostatic tweezer for droplet manipulation. _Proc. Natl. Acad. Sci. U.S.A._ 2022, _119_ , e2105459119
  * Xu et al. 2022 Xu, W.; Jin, Y.; Li, W.; Song, Y.; Gao, S.; Zhang, B.; Wang, L.; Cui, M.; Yan, X.; Wang, Z. Triboelectric wetting for continuous droplet transport. _Sci. Adv._ 2022, _8_ , eade2085
  * Zhang et al. 2022 Zhang, K.; Li, J.; Fang, W.; Lin, C.; Zhao, J.; Li, Z.; Liu, Y.; Chen, S.; Lv, C.; Feng, X.-Q. An energy-conservative many-body dissipative particle dynamics model for thermocapillary drop motion. _Phys. Fluids_ 2022, _34_ , 052011
  * Dundas et al. 2009 Dundas, D.; McEniry, E. J.; Todorov, T. N. Current-driven atomic waterwheels. _Nat. Nanotechnol._ 2009, _4_ , 99–102
  * Regan et al. 2004 Regan, B. C.; Aloni, S.; Ritchie, R. O.; Dahmen, U.; Zettl, A. Carbon nanotubes as nanoscale mass conveyors. _Nature_ 2004, _428_ , 924
  * Zhao et al. 2010 Zhao, J.; Huang, J.-Q.; Wei, F.; Zhu, J. Mass transportation mechanism in electric-biased carbon nanotubes. _Nano Lett._ 2010, _10_ , 4309–4315
  * Kudernac et al. 2011 Kudernac, T.; Ruangsupapichat, N.; Parschau, M.; Maciá, B.; Katsonis, N.; Harutyunyan, S. R.; Ernst, K.-H.; Feringa, B. L. Electrically driven directional motion of a four-wheeled molecule on a metal surface. _Nature_ 2011, _479_ , 208–211
  * Shklyaev et al. 2013 Shklyaev, O. E.; Mockensturm, E.; Crespi, V. H. Theory of carbomorph cycles. _Phys. Rev. Lett._ 2013, _110_ , 156803
  * Fennimore et al. 2003 Fennimore, A. M.; Yuzvinsky, T. D.; Han, W.-Q.; Fuhrer, M. S.; Cumings, J.; Zettl, A. Rotational actuators based on carbon nanotubes. _Nature_ 2003, _424_ , 408–410
  * Bailey et al. 2008 Bailey, S. W. D.; Amanatidis, I.; Lambert, C. J. Carbon nanotube electron windmills: A novel design for nanomotors. _Phys. Rev. Lett._ 2008, _100_ , 256802
  * Huang et al. 2014 Huang, Y.; Zhu, S.; Li, T. Directional transport of molecular mass on graphene by straining. _Extreme Mech. Lett._ 2014, _1_ , 83–89
  * Santos; and Ondarquhus 1995 Santos;, F. D.; Ondarquhus, T. Free-Running Droplets. _Phys. Rev. Lett._ 1995, _75_ , 2972
  * Lee et al. 2002 Lee, S. W.; Kwok, D. Y.; Laibinis, P. E. Chemical influences on adsorption-mediated self-propelled drop movement. _Phys. Rev. E_ 2002, _65_ , 9
  * Daniel and Chaudhury 2002 Daniel, S.; Chaudhury, M. K. Rectified motion of liquid drops on gradient surfaces induced by vibration. 2002, _18_ , 3404–3407
  * Brunet et al. 2007 Brunet, P.; Eggers, J.; Deegan, R. D. Vibration-induced climbing of drops. _Phys. Rev. Lett._ 2007, _99_ , 3–6
  * Brunet et al. 2009 Brunet, P.; Eggers, J.; Deegan, R. D. Motion of a drop driven by substrate vibrations. _Eur. Phys. J.: Spec. Top._ 2009, _166_ , 11–14
  * Kwon et al. 2023 Kwon, O. K.; Kim, J. M.; Kim, H. W.; Kim, K. S.; Kang, J. W. A Study on Nanosensor Based on Graphene Nanoflake Transport on Graphene Nanoribbon Using Edge Vibration. _J. Electr. Eng. Technol._ 2023, _18_ , 663–668
  * Buguin et al. 2002 Buguin, A.; Talini, L.; Silberzan, P. Ratchet-like topological structures for the control of microdrops. _Appl. Phys. A: Mater. Sci. Process._ 2002, _75_ , 207–212
  * Thiele and John 2010 Thiele, U.; John, K. Transport of free surface liquid films and drops by external ratchets and self-ratcheting mechanisms. _Chem. Phys._ 2010, _375_ , 578–586
  * Noblin et al. 2009 Noblin, X.; Kofman, R.; Celestini, F. Ratchetlike motion of a shaken drop. _Phys. Rev. Lett._ 2009, _102_ , 1–4
  * Ni et al. 2022 Ni, E.; Song, L.; Li, Z.; Lu, G.; Jiang, Y.; Li, H. Unidirectional self-actuation transport of a liquid metal nanodroplet in a two-plate confinement microchannel. _Nanoscale Adv._ 2022, _4_ , 2752–2761
  * Barthlott et al. 2016 Barthlott, W.; Mail, M.; Neinhuis, C. Superhydrophobic hierarchically structured surfaces in biology: evolution, structural principles and biomimetic applications. _Phil. Trans. R. Soc. A_ 2016, _374_ , 2016019
  * Badr et al. 2022 Badr, R. G. M.; Hauer, L.; Vollmer, D.; Schmid, F. Cloaking Transition of Droplets on Lubricated Brushes. _J. Phys. Chem. B_ 2022, _126_ , 7047–7058
  * Button et al. 2012 Button, B.; Cai, L.-H.; Ehre, C.; Kesimer, M.; Hill, D. B.; Sheehan, J. K.; Boucher, R. C.; Rubinstein, M. A Periciliary Brush Promotes the Lung Health by Separating the Mucus Layer from Airway Epithelia. _Science_ 2012, _337_ , 937–941
  * Stukowski 2010 Stukowski, A. Visualization and analysis of atomistic simulation data with OVITO–the Open Visualization Tool. _Modelling Simul. Mater. Sci. Eng._ 2010, _18_ , 015012
  * Kremer and Grest 1990 Kremer, K.; Grest, G. S. Dynamics of entangled linear polymer melts: A molecular-dynamics simulation. _J. Chem. Phys._ 1990, _92_ , 5057
  * Theodorakis and Fytas 2011 Theodorakis, P. E.; Fytas, N. G. Microphase separation in linear multiblock copolymers under poor solvent conditions. _Soft Matter_ 2011, _7_ , 1038–1044
  * Tretyakov and Müller 2014 Tretyakov, N.; Müller, M. Directed transport of polymer drops on vibrating superhydrophobic substrates: a molecular dynamics study. _Soft matter_ 2014, _10_ , 4373–86
  * Theodorakis et al. 2010 Theodorakis, P. E.; Paul, W.; Binder, K. Pearl-necklace structures of molecular brushes with rigid backbone under poor solvent conditions: A simulation study. _J. Chem. Phys._ 2010, _133_ , 104901
  * Theodorakis et al. 2010 Theodorakis, P. E.; Paul, W.; Binder, K. Interplay between Chain Collapse and Microphase Separation in Bottle-Brush Polymers with Two Types of Side Chains. _Macromolecules_ 2010, _43_ , 5137–5148
  * Schneider and Stoll 1978 Schneider, T.; Stoll, E. Molecular-dynamics study of a three-dimensional one-component model for distortive phase transitions. _Phys. Rev. B_ 1978, _17_ , 1302–1322
  * Plimpton 1995 Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. _J. Comp. Phys._ 1995, _117_ , 1–19
  * Barber et al. 1996 Barber, C. B.; Dobkin, D. P.; Huhdanpaa, H. T. The Quickhull algorithm for convex hulls. _ACM Trans. on Mathematical Software_ 1996, _22_ , 469–483, http://www.qhull.org
  * Theodorakis et al. 2015 Theodorakis, P. E.; Müller, E. A.; Craster, R. V.; Matar, O. K. Modelling the superspreading of surfactant-laden droplets with computer simulation. _Soft Matter_ 2015, _11_ , 9254–9261