# The Fractal Geometrical Properties of Nuclei

W. H. Ma<EMAIL_ADDRESS>Institute of Modern Physics, Chinese Academy of
Science, Lanzhou 730000, People’s Reublic of China  University of Chinese
Academy of Science, Beijing, 100049, People’s Reublic of China  J. S. Wang
Institute of Modern Physics, Chinese Academy of Science, Lanzhou 730000,
People’s Reublic of China  Q. Wang  Institute of Modern Physics, Chinese
Academy of Science, Lanzhou 730000, People’s Reublic of China  S. Mukherjee
Physics Department, Faculty of Science, M.S. University of Baroda, Vadodara -
390002, India  L. Yang  Institute of Modern Physics, Chinese Academy of
Science, Lanzhou 730000, People’s Reublic of China  Y. Y. Yang  Institute of
Modern Physics, Chinese Academy of Science, Lanzhou 730000, People’s Reublic
of China  M. R. Huang  Institute of Modern Physics, Chinese Academy of
Science, Lanzhou 730000, People’s Reublic of China  Y. J. Zhou  Institute of
Modern Physics, Chinese Academy of Science, Lanzhou 730000, People’s Reublic
of China  University of Chinese Academy of Science, Beijing, 100049, People’s
Reublic of China

###### Abstract

We present a new idea to understand the structure of nuclei, which is
comparing to the liquid drop model. After discussing the probability that the
nuclear system may be a fractal object with the characteristic of self-
similarity, the nuclear irregular structure properties and the self-similarity
characteristic are considered to be an intrinsic aspects of nuclear structure
properties. For the description of nuclear geometric properties, nuclear
fractal dimension is an irreplaceable variable similar to the nuclear radius.
In order to determine these two variables, a new nuclear potential energy
formula which is related to the fractal dimension is put forward and the
phenomenological semi-empirical Bethe-Weizs$\ddot{a}$cker binding energy
formula is modified using the fractal geometric theory. And one important
equation set with two equations is obtained, which is related to the
conception that the fractal dimension should be a dynamical parameter in the
process of nuclear synthesis. The fractal dimensions of the light nuclei are
calculated and their physical meanings are discussed. We have compared the
nuclear fractal mean density radii with the radii calculated by the liquid
drop model for the light stable and unstable nuclei using rational nuclear
fractal structure types. In the present model of fractal nuclear structure
there is an obvious feature comparing to the liquid drop model, since the
present model can reflect the geometric informations of the nuclear structure,
especially for the nuclei with clusters, such as the $\alpha$-cluster nuclei
and halo nuclei.

###### pacs:

##  I Introduction

Both in theory bibitem1 ; bibitem2 ; bibitem3 ; bibitem4 ; bibitem5  and
experiment bibitem6 ; bibitem7 ; bibitem8 ; bibitem9 ; bibitem10 ; bibitem11 ;
bibitem12 , the concept of the cluster structure plays an important role in
nuclear reaction, nuclear structure and the nuclear excitation. The nuclear
distribution is non-uniform especially for halo nuclei. For instance,
${}^{11}Li$ has a core and valence nucleons (${}^{9}Li+n+n$) halo structure
and the matter distribution is much non-uniform bibitem10 ; bibitem13 . Its
radius is much larger than given by the usual expression
($r_{0}A^{\frac{1}{3}}$), which depends on the uniform-density liquid drop
model.

In addition, the phenomenological semi-empirical Bethe-Weizs$\ddot{a}$cker
binding energy formula bibitem14 ; bibitem15 ; bibitem16  for the masses of
nuclei has been derived from the liquid drop model and successfully used to
calculate the binding energy for stable nuclei and the nuclei very near the
stable line. However, there are two points we cannot ignore. On one hand, the
liquid drop model for calculating nuclear binding energy strictly depends on
the experimental data. On the other hand, it is difficult to describe the
light and halo nuclei because it depends on statistic and cannot reflect the
structure properties in nucleus precisely. Taking into consideration the fact
that, due to the nuclear particle property and its quantum motion, it is not
proper to regard nucleus as a compact sphere with smooth surface and there is
no explicit border. The nucleons in a nucleus are separated and there are void
spaces among them. These void spaces near the surface are much bigger in halo
or weakly-bounded nuclei in comparison to that in stable nuclei. The notion of
radii is only a statistical average effect of nuclear matter distribution. In
fact, the nuclear structures are irregular, which are the obvious features
different from the description of the liquid drop model.

In order to describe the nuclear structure properties reasonably well, the
notion of fractal object for the irregular systems with non-integral
dimensions is possible, which was introduced by Mandelbrot in science in 1967
bibitem17 . The self-similarity and the scale invariance are the most
important characteristics of a fractal object. The fractal self-similarity
means that the part of the fractal object is similar to the whole fractal
object after magnification. For an arbitrary part of the fractal object, its
some important properties, such as shape, complexity and irregularity, keep
invariant after magnifying or shrinking certain proportion, which is the
meaning of scale invariance. For a regular fractal object, there is no
characteristic size but characteristic fractal dimension. However, the self-
similarity of an approximate or statistical fractal object exists in a finite
scale range. In nuclear system with clusters structure, due to the similar
nuclear and electromagnetic interaction between the clusters in a nucleus, it
is possible that the geometric structure and the physical laws are similar
between the nucleons in cluster and the clusters in nucleus and also between
the clusters and the whole nucleus. The nucleus can be deemed as a statistical
self-similarity fractal system with finite scale range. Xavier Campi had found
the existence of finite size scaling in nuclear fragmentation bibitem18 . He
introduced the fractal dimension which give the information on the internal
structure of the fragments during the fragmentation process. Adamenko et al
bibitem19  queried the liquid-like model and studied the properties in super-
heavy nuclear isomers using fractal theory. These studies are the fine
enlightenment of connecting the characteristics of the nuclear structure with
the fractal geometric theory.

As explained above, the nuclear irregular structure properties and the self-
similarity characteristic may be an intrinsic aspects of nuclear structure
properties. It is possible to use the fractal theory to describe the nuclear
structure properties. The objective of the present work is to introduce this
new concept of nuclear fractal structures to study the nuclear properties.

##  II Description of the model of nuclear fractal structure

In our present description, a more general conception of nuclear fractal
clusters (NFCs) is applied, which is similar but different from the
conventional one of the $\alpha$-cluster structure and the core plus valence
nucleons structure in halo nuclei. The latter is considered to be one kind of
the former which is basing on the concept of the characteristics of fractal
objects. The concept of NFCs is that of the nucleus as a fractal assembly of
structural subunits that are themselves made up of no less than one nucleon
and keep certain correlation of the similarity with the ensemble in the
geometrical and in physics. And the geometrical boundary of the NFCs within
some nuclei is less distinguishable than that of the clusters in
$\alpha$-cluster nuclei and the halo nuclei.

Similar to the definition given in bibitem20 , an isotropic self-similarity
nuclear fractal dimension $D_{f}$ is defined in the following relation

$M(b\cdot{r})=b^{D_{f}}\cdot{M(r)}$ (1)

where $M(r)$ is the mass number within the size $r$ of the fractal object;
$M(b\cdot{r})$ is the mass number of $b$ times the size $r$ of the fractal
object, where $b$ is a scaling factor among the similar parts within the
fractal object. The only solution for relation (1) is
$M(r)\propto{r^{D_{f}}}$. The nuclear average matter density $\rho(r)$ with
the law of decay of isotropic spatial correlation,
$\rho(r)=\frac{M(r)}{V(r)}\varpropto{\frac{r^{D_{f}}}{r^{3}}}\varpropto{r^{D_{f}-3}}$
, is a basic variable function in nucleus. So far, the geometric dimension of
nuclei is considered as 3, because of the concept of liquid drop model. For a
real physical nuclear object embedded in 3-dimension Euclidean-space, its
dimension must be less than or equal to 3. Most of the nuclei with fractal
dimensions approaching 3 are stable, which are more like liquid drops.

In relation (2), we assume that there are several NFCs and $A_{i}$ is the mass
number of a NFC within a nucleus with mass number $A$. $\rho_{dis}(r)$ is the
density distribution function. $\rho_{i}$ is the fractal mean density of a NFC
and $R_{i}$ is the fractal mean density radius of it. $F$ is the number of the
NFCs in a given nucleus.

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\frac{A}{V}=\frac{\int{\rho_{dis}(r)\cdot{dv}}}{\int{dv}}=\frac{\sum_{i=1}^{F}A_{i}}{V}$
(2) $\displaystyle=$
$\displaystyle\frac{\sum_{i=1}^{F}\int\rho_{dis}(r_{i})\cdot{dv_{i}}}{V}=\frac{\sum_{i=1}^{F}\rho_{i}\cdot{V_{i}}}{V}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{F}\frac{V_{i}}{V}\cdot{\rho_{i}}=\sum_{i=1}^{F}(\frac{R_{i}}{R})^{3}\cdot{\rho_{i}}$

where, $R$ is the nuclear fractal mean density radius defined by the relation
$A=\frac{4}{3}\pi{R^{3}}\rho$. And $\rho_{i}=\frac{3A_{i}}{4\pi{R_{i}^{3}}}$.

In addition, due to the assumption of nuclear fractal structure and the mass-
radius relation of self-similarity fractal object, we list several basic
relations:

$A=\sum_{i=1}^{F}A_{i};$ (3)
$A\propto{R^{D_{f}}};A_{i}\propto{R_{i}^{D_{f}}}.$ (4)

Using relations (3) and (4), we get the relation among $R$, $R_{i}$ and
$D_{f}$:

$R^{D_{f}}=\sum_{i=1}^{F}R_{i}^{D_{f}}.$ (5)

Due to the two proportional relations in (4), the relation among $R$, $R_{i}$,
$A$, $A_{i}$ and $D_{f}$ is

$R=R_{i}(\frac{A}{A_{i}})^{\frac{1}{D_{f}}}.$ (6)

Substitute (6) into (2), then get

$\rho=\sum_{i=1}^{F}(\frac{A_{i}}{A})^{\frac{3}{D_{f}}}\rho_{i}.$ (7)

And the relation between the NFCs in nucleus is

$\rho_{i}=\rho_{j}(\frac{A_{i}}{A_{j}})^{\frac{D_{f}-3}{D_{f}}}$ (8)

Therefore, to describe the nuclear geometric properties, the nuclear fractal
dimension is an irreplaceable variable similar to the nuclear radius. For
given nuclear system with cluster structure, the final free variables are
($R$, $D_{f}$) or ($\rho$, $D_{f}$ ), which are also ($R_{i}$, $D_{f}$) or
($\rho_{i}$, $D_{f}$) due to the relations (6) and (7). In order to determine
these two variables, we put forward a new nuclear potential energy formula
which is related to fractal dimension. Then the phenomenological semi-
empirical Bethe-Weizs$\ddot{a}$cker binding energy formula is modified and the
total potential energy is obtained. And one important equation set with two
equations is obtained, which is related to the conception that the fractal
dimension should be a dynamical parameter in the process of nuclear synthesis.
So, the calculations of nuclear fractal dimensions and radii can be done.

In nuclear system, considering the self-similarity properties of the nuclear
fractal system, we put forward a nuclear potential energy formula

$u(r)=\frac{v_{0}D_{f}}{3(D_{f}-2)}(\frac{r}{r_{s}})^{D_{f}-3},2<D_{f}\leq{3}.$
(9)

Which is proportional to the nuclear average density $\rho(r)$. $2r_{s}$stands
for the minimum scale size of a nuclear fractal system and it is also the
maximum size of the minimum cluster element. $v_{0}$ is a coefficient and
keeps constant, which is also corresponding to the estimation of the depth of
nuclear potential well in liquid drop model when $D_{f}=3$.

From (9) and using the idea of self-similarity, the relation for the nuclear
potential energy of one NFC within nucleus is obtained:

$u_{i}(r_{i})=\frac{v_{0}D_{f}}{3(D_{f}-2)}(\frac{r_{i}}{r_{s}})^{D_{f}-3},2<D_{f}\leq{3}$
(10)

And because of $u(R)\propto{\rho(R)}$ ,$u_{i}(R_{i})\propto{\rho_{i}(R_{i})}$
and (7), the relation between $u(R)$ and $u_{i}(R_{i})$ is

$u(R)=\sum_{i=1}^{F}(\frac{A_{i}}{A})^{\frac{3}{D_{f}}}u_{i}(R_{i}).$ (11)

Next, we modify the phenomenological semi-empirical Bethe-Weizs$\ddot{a}$cker
binding energy formula with the fermi gas model and the fractal theory. Here
we mainly concern the liquid drop energy and put aside the correction term
basing on the microscopic method, such as the description in bibitem21 . The
original one derived from liquid drop model is

$\displaystyle
B=(u_{depth}-c_{v}-c_{as}(1-\frac{2Z}{A})^{2})A-c_{surf}A^{\frac{2}{3}}$ (12)
$\displaystyle-
c_{Q}\frac{Z(Z-1)}{A^{\frac{1}{3}}}+c_{p}\frac{(-1)^{Z}+(-1)^{A-Z}}{2A^{\frac{4}{3}}};$

where $u_{depth}\approx{58}$ bibitem22 (the estimation of the depth of nuclear
potential well); $c_{v}=42.27$; $c_{as}=23.48$; $c_{surf}=17.72$;
$c_{Q}=0.72$; $c_{p}=19.39$. We use the experimental mass data bibitem23  to
fit the other parameters.

The modified formulae are:

$B_{strong}=(v_{depth}-c_{1}(\rho)-c_{2}(\rho)(1-\frac{2Z}{A})^{2})A;$ (13)
$B_{surf}=-c_{s}4\pi{R^{2}};$ (14) $B_{Q}=-\frac{3}{5}\frac{Z(Z-1)e^{2}}{R};$
(15) $B_{p}=c_{p}\frac{(-1)^{Z}+(-1)^{A-Z}}{2A^{\frac{4}{3}}};$ (16)
$B=B_{strong}+B_{surf}+B_{Q}+B_{p};$ (17)

where, $c_{1}(\rho)=\frac{3}{5}\varepsilon(\rho)\rho^{\frac{2}{3}}$;
$c_{2}(\rho)=\frac{1}{3}\varepsilon(\rho)\rho^{\frac{2}{3}}$;
$c_{s}=0.98MeV\cdot{fm^{2}}$; $c_{p}=19.39MeV$;
$\varepsilon(\rho)=\varepsilon_{0}(\frac{\rho}{\rho_{0}})^{\alpha}$;
$\varepsilon_{0}=264.12MeV^{-1}$; $\rho_{0}=0.138fm^{-3}$. The formula (13) is
the volume energy consisting of three parts. The first part is the depth of
nuclear potential energy, which is not a constant, but depends on the
structure type. The second part and the third part derived from fermi-gas
model are the kinetic energy of the nuclear system, because of which the
binding energy decrease. $\varepsilon(\rho)$ shares the kinetic energy, which
decrease the contribution of kinetic energy. $\varepsilon_{0}$ was explained
as a constant related to the virial coefficient in bibitem19 . Since we
consider that the nuclear system has fractal structure that is much different
from the fermi gas, $\varepsilon_{0}$ cannot be kept constant any more, but
depends on the mean density that is related to the mean density of every
cluster in nucleus as relation (7) shows. In the formula of
$\varepsilon(\rho)$, $\alpha$ is a parameter. The formulas(14) and (15) for
surface energy and the energy because of Coulomb interaction dose not have any
parameters. $c_{s}$, the coefficient of surface tension, is about
$1MeV/fm^{2}$. Here we still keep the pair energy because the binding energy
is larger when the proton and neutron numbers are even and it is smaller when
one of the numbers are odd and even more so when both are odd, which maybe has
nothing to do with the nuclear fractal structure. Thus all the terms
contributing to binding energy are depended on the fractal dimension $D_{f}$
except the pair term. When $D_{f}\rightarrow{3}$ , the modified formula
degenerates to the ordinary binding energy formula (12), and then the state of
nucleus shifts from the non-uniform density, which have NFCs, to the uniform
state of liquid drop model. So, the parameters $c_{s}$, $c_{p}$,
$\varepsilon_{0}$ and $\rho_{0}$ can be all derived from (12-17) When
$D_{f}\rightarrow{3}$.

A possible formula for the depth of nuclear potential well in the nuclear
fractal system is found to be

$\displaystyle v_{depth}=\frac{1}{2}(1-(\frac{A-(1+s)Z}{A})^{2})$
$\displaystyle\times{\sum_{i=1}^{F}\frac{A_{i}}{A}(58+3u_{i}(R_{i})\frac{D_{f}-2}{D_{f}})}$
(18)

which depends on the nuclear fractal structure. In fact, the depth of nuclear
potential well cannot be $58Mev$ for all nuclei because of the different
fractal structure for different nuclei. The product factor is introduced to
describe the difference between the number of neutrons and protons in nucleus,
and $s=\frac{N_{s}}{Z_{s}}$, $N_{s}$ and $Z_{s}$ are the number of neutrons
and protons in nucleus on stable line, which makes sure the depth of nuclear
potential well get the maximum value for the nuclei on stable line. For the
light stable nuclei, $s\approx{1}$. And for the nuclei ${}^{4}He$, ${}^{8}Be$,
${}^{12}C$, ${}^{16}O$, ${}^{20}Ne$, ${}^{24}Mg$, ${}^{28}Si$, ${}^{32}S$,
${}^{36}Ar$ and ${}^{40}Ca$, $s=1$. When $D_{f}\rightarrow{3}$,
$u=u_{i}=v_{0}=-u_{depth}.$

Figure 1:  The interaction among clusters in given nucleus. The notations i,
j, k stands for different NFCs.

We define the total potential energy $U=U(A,Z,D_{f},\rho)$, which is sum of
the total nuclear potential energy $U_{A}=U_{A}(A,D_{f},\rho)$ and the total
electromagnetic potential energy $U_{Z}=U_{Z}(Z,D_{f},\rho)$. Namely,

$U=U_{A}+U_{Z}.$ (19)

The interaction among clusters in given nucleus is showed in FIG. 1. The total
nuclear potential energy is

$U_{A}=\sum_{i=1}^{F}(\frac{A_{i}}{A})^{\frac{3}{D_{f}}}U_{i}+\sum_{i=1}^{F}\sum_{i\neq{j}}^{F}(\frac{A_{j}}{A})^{\frac{3}{D_{f}}}U_{ij},$
(20)

Where,

$\displaystyle U_{i}$ $\displaystyle=$ $\displaystyle
4\pi\frac{v_{0}D_{f}}{3(D_{f}-2)}\int_{0}^{R_{i}}(\frac{r_{i}}{r_{s}})^{D_{f}-3}\rho_{dis}(r_{i})r_{i}^{2}dr_{i}$
(21) $\displaystyle=$
$\displaystyle\frac{4\pi{v_{0}\rho_{i}R_{i}^{D_{f}}}}{3(D_{f}-2)r_{s}^{D_{f}-3}}.$

$\displaystyle U_{ij}$ $\displaystyle=$ $\displaystyle
2{\pi}\int_{0}^{\pi}\int_{0}^{R_{j}}u_{i}(r)\rho_{dis}(r_{j})r_{j}^{2}\sin(\theta)dr_{j}d\theta$
$\displaystyle=$ $\displaystyle
2\pi\frac{v_{0}D_{f}}{3(D_{f}-2)}\frac{\rho_{j}(R_{ij}+R_{j})^{D_{f}+1}G(D_{f})}{(D_{f}+D_{f}^{3})R_{ij}(R_{ij}^{2}-R_{j}^{2})r_{s}^{D_{f}-3}}.$

In(22),
$G(D_{f})=(R_{ij}^{2}+D_{f}R_{j}^{2})((\frac{R_{ij}-R_{j}}{R_{ij}+R_{j}})^{D_{f}+1}-1)+(D_{f}+1)R_{ij}R_{j}((\frac{R_{ij}-R_{j}}{R_{ij}+R_{j}})^{D_{f}+1}+1)$;
the distance $r=\sqrt{r_{j}^{2}+R_{ij}^{2}-2r_{j}R_{ij}\cos{\theta}}$ as
showed in FIG. 1; $R_{ij}\approx{R_{i}+R_{j}}.$ The total electromagnetic
potential energy is

$U_{Z}=\sum_{i=1}^{F}U_{Z_{i}}+\sum_{i=1}^{F}\sum_{i\neq{j}}^{F}U_{Z_{ij}}$
(23)

where,

$U_{Z_{i}}=\frac{3}{5}Z_{i}(Z_{i}-1)\frac{e^{2}}{R_{i}};$ (24)
$U_{Z_{ij}}=\frac{Z_{i}Z_{j}e^{2}}{R_{ij}}$ (25)

For a given nucleus with NFCs in it, we consider that its fractal dimension
$D_{f}$ is constant. However it should be a dynamical parameter in the process
of nuclear synthesis. In a given nuclear reaction the binding energy
$B=B(A,Z,D_{f},\rho)$ and the total potential energy $U=U(A,Z,D_{f},\rho)$ are
changing with $D_{f}$. When $D_{f}$ gets the fixing value, then
$B=B(A,Z,D_{f},\rho)$ and $U=U(A,Z,D_{f},\rho)$ get the minimum value, which
corresponds to an interacting system becoming relatively stable. So one
equation set is gotten:

$\left\\{\begin{array}[]{c}\partial_{D_{f}}B(A,Z,D_{f},\rho)=0\\\
\partial_{D_{f}}U(A,Z,D_{f},\rho)=0.\end{array}\right.$ (26)

Because of the uncertain parameter $\alpha$, we need more than two equations
to study the structure properties in nucleus. The additional one is

$B(A,Z,D_{f},\rho)=B_{exp}+E_{excited},$ (27)

where $B_{exp}$ is the experimental value of binding energy. And the
$E_{excited}$ is the change of the binding energy due to the nucleus being
excited from the ground state, which is corresponding to the situation that
the NFCs structure is forming in the excited nucleus. If only the nuclei in
ground states are considered, $E_{excited}=0$.

Finally, we arrive at the modified binding energy formula $B(A,Z,D_{f},\rho)$
and the total potential energy $U(A,Z,D_{f},\rho)$, which are the functions of
$A$, $Z$, $D_{f}$ and $\rho$. And the important equation set (26) is gotten.
Use (26), (27) and assume rational fractal structure types in nucleus, the
nuclear fractal dimension and radius can be determined.

##  III Results and discussions

In the present work, we have done some calculations and obtained some results
for light nuclei in ground states. Combining the analyses given in bibitem1 ;
bibitem2 ; bibitem3 ; bibitem4 ; bibitem5 ; bibitem6 ; bibitem7 ; bibitem8 ;
bibitem9 ; bibitem10 ; bibitem11 ; bibitem12  and bibitem24 ; bibitem25 ;
bibitem26  and the consideration that the rational fractal structure of one
nucleus is determined by making the seperated energy of NFCs be minimum
possibly, the possible NFC structures for these nuclei are shown in TABLE I.
We consider that the scale $r_{s}$ is about 1.25 fm, which is the half of the
maximum size of one nucleon. And we compare the radii
($R_{L}=r_{0}A^{\frac{1}{3}}$,$r_{0}=1.25fm$) from liquid drop model with the
nuclear fractal mean density radii.

Nuclei | NFCs | Nuclei | NFCs
---|---|---|---
${}^{5}He$ |  ${}^{4}He$+n | ${}^{14}B$ |  ${}^{13}B$+n
${}^{6}He$ |  ${}^{4}He$+n+n | ${}^{15}B$ |  ${}^{14}B$+n
${}^{7}He$ |  ${}^{6}He$+n | ${}^{16}B$ |  ${}^{15}B$+n
${}^{8}He$ |  ${}^{6}He$+n+n | ${}^{9}C$ |  ${}^{8}B$+p
${}^{5}Li$ |  ${}^{4}He$+p | ${}^{10}C$ |  ${}^{8}Be$+p+p
${}^{6}Li$ |  ${}^{4}He$+${}^{2}H$ | ${}^{11}C$ |  ${}^{7}Be$+${}^{4}He$
${}^{7}Li$ |  ${}^{6}Li$+n | ${}^{12}C$ |  ${}^{8}Be$+${}^{4}He$
${}^{8}Li$ |  ${}^{7}Li$+n | ${}^{13}C$ |  ${}^{12}C$+n
${}^{9}Li$ |  ${}^{8}Li$+n | ${}^{14}C$ |  ${}^{13}C$+n
${}^{10}Li$ |  ${}^{9}Li$+n | ${}^{15}C$ |  ${}^{14}C$+n
${}^{11}Li$ |  ${}^{9}Li$+n+n | ${}^{16}C$ |  ${}^{15}C$+n
${}^{8}Be$ |  ${}^{4}He$+${}^{4}He$ | ${}^{13}N$ |  ${}^{12}C$+p
${}^{9}Be$ |  ${}^{8}Be$+n | ${}^{14}N$ |  ${}^{10}B$+${}^{4}He$
${}^{10}Be$ |  ${}^{9}Be$+n | ${}^{15}N$ |  ${}^{11}B$+${}^{4}He$
${}^{11}Be$ |  ${}^{10}Be$+n | ${}^{16}N$ |  ${}^{15}N$+n
${}^{12}Be$ |  ${}^{11}Be$+n | ${}^{17}N$ |  ${}^{16}N$+n
${}^{13}Be$ |  ${}^{12}Be$+n | ${}^{13}O$ |  ${}^{12}N$+p
${}^{14}Be$ |  ${}^{12}Be$+n+n | ${}^{14}O$ |  ${}^{13}N$+p
${}^{8}B$ |  ${}^{7}Be$+n | ${}^{15}O$ |  ${}^{14}N$+p
${}^{9}B$ |  ${}^{8}Be$+n | ${}^{16}O$ |  ${}^{12}C$+${}^{4}He$
${}^{10}B$ |  ${}^{6}Li$+${}^{4}He$ | ${}^{17}O$ |  ${}^{13}C$+${}^{4}He$
${}^{11}B$ |  ${}^{7}Li$+${}^{4}He$ | ${}^{18}O$ |  ${}^{14}C$+${}^{4}He$
${}^{12}B$ |  ${}^{11}B$+n |  |
Table 1: The possible NFC structures for light nuclei. Nuclei | $R_{f}(fm)$ | $R_{m}(fm)$ | Ref.
---|---|---|---
${}^{6}Li$ | $2.217$ | $2.32\pm{0.03}$ | bibitem27
${}^{6}He$ | $2.583$ | 2.45(10) | bibitem28
${}^{7}Li$ | $2.374$ | $2.33\pm{0.02}$ | bibitem27
${}^{8}He$ | $3.100$ | $2.53(8)$ | bibitem28
${}^{8}B$ | $2.635$ | $2.55\pm{0.08}$ | bibitem29
${}^{8}Li$ | $2.629$ | $2.583\pm{0.023}$ | bibitem27
${}^{9}C$ | $2.874$ | $2.71(32)$ | bibitem30
${}^{10}B$ | $2.647$ | $2.56\pm{0.23}$ | bibitem27
${}^{11}Li$ | $3.432$ | $3.34_{-0.08}^{+0.04}$ | bibitem31
${}^{11}Be$ | $3.027$ | $3.039\pm{0.038}$ | bibitem32
${}^{12}B$ | $3.016$ | $2.723\pm{0.050}$ | bibitem27
${}^{12}C$ | $2.848$ | $2.48\pm{0.08}$ | bibitem27
${}^{13}B$ | $3.193$ | $2.746\pm{0.050}$ | bibitem27
${}^{14}Be$ | $3.736$ | $3.36\pm{0.19}$ | bibitem27
${}^{14}B$ | $3.383$ | $3.00\pm{0.10}$ | bibitem27
${}^{14}N$ | $3.027$ | $2.61\pm{0.10}$ | bibitem27
${}^{16}O$ | $3.194$ | $2.63\pm{0.06}$ | bibitem27
Table 2: The comparison of the calculated values of the nuclear fractal mean
density radii ( here we use $R_{f}$ to stand for the nuclear fractal mean
density radii for convenience.) using the model of nuclear fractal structure
with the experimental results of root-mean-square (rms) matter ($R_{m}$)
radii. The cluster structures of the nuclei listed in here are same as those
in TABLE I.

In order to describe the essential feature of the nuclear fractal self-
similarity symmetry, fractal dimension is one of the basic geometrical
parameters. It s related to the nuclear homogeneity and the ingredients in
nucleus, which reflect the nuclear irregular degree. It can be seen that the
values of the dimensions are much close to 3, if the mean densities of every
NFC in one nucleus and the mean density of this nucleus are all not much
different and the irregular degrees of the nuclear structures are low. Based
on the NFC structures (TABLE I) and the values of the dimensions (FIG. 2),
some interesting results may be discussed. For isobars, because the degree of
homogeneity of the stable nuclei is higher than the nuclei far from the stable
line generally, the values of the fractal dimension for the former are greater
than for the latter. The fractal dimensions of isobars are comparative or
same, if the NFC structure types of these nuclei are similar. For example, the
fractal dimensions of the isobars with two-body structure type, (${}^{5}Li$,
${}^{5}He$), (${}^{8}Li$, ${}^{8}Be$, ${}^{8}B$), (${}^{10}Be$, ${}^{10}B$),
(${}^{11}C$, ${}^{11}B$), (${}^{13}N$, ${}^{13}C$) and (${}^{16}C$,
${}^{16}N$), are almost same respectively. For nuclei with same $A$, the
fractal dimensions of the nuclei with three-body structure, such as,
${}^{6}He$, ${}^{8}He$, ${}^{10}C$, ${}^{11}Li$ and ${}^{14}Be$, are
distinctly lower than the nuclei with two body structure type. The values of
fractal dimension decrease due to the increase of nuclear mass number, which
may be associated with the scale $r_{s}$. The fractal scale variables amount
to the resolution of the measurement bibitem33 , which are related to the
scale relativity. For instance, if the half of the maximum size of one nucleon
serves as the fractal scale, using a 3-dimension bulk whose radius equals to
this scale to cover a nucleon, the dimension of one nucleon is 3. However,
along with the increase of the nuclear mass number, the nuclear structure
becomes more irregular. Therefore, using this 3D structure to measure the
nucleus, the dimension of this nucleus is lower than 3, which reflects the
information of the nuclear irregular structures. In brief, as the results of
the effects of all structure variables and binding energy, the fractal
dimension can generally describe the nuclear structure well. Besides, in the
process of nuclear synthesis $D_{f}$ is a dynamical parameter, which when
$D_{f}$ gets a fixed value, then it explains the existence of relatively
stable nuclei.

Figure 2:  the values of the fractal dimensions of the nuclei. Figure 3:  the
comparison of the calculated values of radii using the model of nuclear
fractal structure with the liquid drop model. The color symbols stand for the
nuclear fractal mean density radii.

The nuclear fractal mean density radii depending on the scale variable $r_{s}$
in the process of measurement can represents the nuclear size statistically,
which are mainly determined by the binding energy and NFC structure types and
related to the fractal dimension by the expression (5). For isobars, the more
is the binding energy, the less is the radii. Generally, as showed in FIG. 3,
for stable nuclei and the nuclei very near the stable line, the values of the
nuclear fractal mean density radii are approaching the line of $R_{L}$
depending on liquid drop model. When the nuclei are far from the line of
stability, the values of the nuclear fractal mean density radii will be
greater than $R_{L}$ and this will increase as the deviation from the line of
stability increases. In some special cases, the nuclei having same NFC
structure and having almost same binding energy, such as ${}^{5}Li$ and
${}^{5}He$, the radii are nearly same. The same situations arises in the pair
of isobars (${}^{8}B$, ${}^{8}Li$), (${}^{9}Be$, ${}^{9}B$) and (${}^{13}C$,
${}^{13}N$). As a matter of fact, the liquid drop model cannot calculate the
radii of the nuclei far from stable line successfully, as they have much
larger radii contrary of the predictions by liquid drop model. On the
contrary, the present fractal structure model can successfully predict the
radii of such nuclei far from the line of stability. Moreover, the comparison
of the calculated values of the nuclear fractal mean density radii using the
model of nuclear fractal structure with the experimental results of root-mean-
square (rms) matter ($R_{m}$) radii are listed in TABLE II, where,
$R_{m}=\sqrt{\frac{Z}{A}R_{p}^{2}+\frac{N}{A}R_{n}^{2}}$ bibitem31 . And, for
the stable nuclei, the results of the calculated radii both using the liquid
drop model and the present method are larger than the experimental results of
root-mean-square (rms) matter ($R_{m}$).

With regard to the parameter $\alpha$ introduced in the second part, it can be
solved as same time as the fractal dimension $D_{f}$ and the nuclear fractal
mean density radius $R$ through the three equations explained in the end of
the second part, which vary from 0 to 1 in present calculations generally
corresponding to the nuclei from stability to unstability.

In summary, the NFC structures are determined by the interactions within
nuclear systems and that the fractal dimension can generally describe such
structural features. The nuclear fractal mean density radii represent
approximately the nuclear size, which are associated with the scale variables.
Actually, the relations among the nuclear structure geometric variables are
complex and correlative.

##  IV Conclusions

In the present work, we consider the importance that the nuclear irregular
structure properties and the self-similarity characteristic may be the
intrinsic aspects of nuclear structure properties. For the description of
nuclear geometric properties, nuclear fractal dimension is an irreplaceable
variable similar to the nuclear radius. Compared with the liquid drop model,
it is a feature that the fractal description can reflect the important
characteristics of the NFC structures especially for describing the nuclei far
from the line of stability and $\alpha$-cluster nuclei. Similar to the liquid
drop model, the present model can get the same results of the relation between
$A$ and $Z$ on the $\beta$ stable line for the light stable nuclei and cannot
predict the existence of the magic nuclei. For heavier nuclei, more realistic
formulae for the depth of nuclear potential energy are needed to be obtained
whose fractal structure is more complicated. These heavy nuclei may have
multilevel fractal structure with several clusters in every level. An
anisotropic description for the nuclear fractal structure and the scale
dependent properties of the nuclear fractal system would be further studied in
future. Further studies will be focused on combination between the properties
of the nuclear fractal structure and the quantum mechanics.

##  V Acknowledgments

This work was partially supported by the Major State Basic Research
Development Program of China (973 Program: New physics and technology at the
limits of nuclear stability) and the Directed Program of Innovation Project of
the Chinese Academy of Sciences with Grant No. KJCX2-YW-N44. The author
Wei.Hu.Ma thanks the University of Chinese Academy of Sciences for the
education. And one of the authors S.Mukherjee acknowledges the financial help
from TWAS-UNESCO, Italy, for visiting Institute of Modern Physics-Chinese
Academy of Science, Lanzhou, China for this work.

## References

## References

  * (1) J. -P. Ebran, et al., Nature. 487, 341 (2012).
  * (2) T. Tomoda, A. Arima, Nucl. Phys. A 303, 217 (1978).
  * (3) M. K. Basu, Phys. Rev. C 6, 2 (1972).
  * (4) R. Beck, R. Krivec, M.V. Mihailovi$\acute{c}$, Nucl. Phys. A 363, 365 (1981).
  * (5) R. A. Baldock B.Buck, J. Phys. G 12, 1 (1986).
  * (6) T. A.Tombrello, G. C.Phillips, Phys. Rev. 122, 1 (1961).
  * (7) P.Taras, et al., Phys. Lett. B 107, 1 (1981).
  * (8) T.Yamaya, et al., Phys. Rev. C 42, 5 (1990).
  * (9) B. Buck, A. C. Merchant, S. M. Perez, Phys. Rev. Lett. 72, 9 (1994).
  * (10) M. V. Zhukov, et al., Phys. Rev. C 44, 1 (1991).
  * (11) L. Giot, et al., Phys. Rev. C 71, 064311 (2005).
  * (12) R. Palit, et al., Phys. Rev. C 68, 034318 (2003).
  * (13) I. Tanihata et al., Phys. Rev. Lett 55, 24 (1985).
  * (14) C. F. von Weizs$\ddot{a}$cker, Z. Physik 96, 461 (1935).
  * (15) H. A. Beth, F. Bacher, Rev. Mod. Phys. 8, 82 (1936).
  * (16) G. Royer, Nucl. Phys. A 807, 105 (2008).
  * (17) B. Mandelbrot, Science. 156, 3775 (1967).
  * (18) X. Campi, Nucl.Phys. A 495, 259 (1989).
  * (19) S. Adamenko, V. Bolotov, V. Novikov, Interdisciplinary Studies of Complex Systems. 1, 1 (2012).
  * (20) Benoit B. Mandelbrot, Phys. Scr. 32, 257 (1985).
  * (21) N. Wang and M. Liu, J.Phys.Conf.Series 420, 012057 (2013).
  * (22) A. B. Migdal, et al., Phys. Rep. 192, 4 (1990).
  * (23) M. Wang, et al., Chinese.Phys. C 36, 12 (2012).
  * (24) W. von Oertzen, et al., Phys. Rep. 432, 43 (2006).
  * (25) A. S. Jensen, et al., Rev. Mod. Phys. 76, 1 (2004).
  * (26) A. H. Wuosmaa, et al., Annu. Rev. Nucl. Part. Sci. 45, 89 (1995).
  * (27) A. Ozawa, et al., Nucl. Phys. A 693, 32 (2001).
  * (28) G. D. Alkhuxx, et al., Nucl. Phys. A 734, 361-364 (2004).
  * (29) S. S. Chandel and S. K. Dhiman, Phys. Rev. C 68, 054320 (2003).
  * (30) R. E. Warner, et al., Phys. Rev. C 74, 014605 (2006).
  * (31) T. Moriguchi, et al., Phys. Rev. C 88, 024610 (2013).
  * (32) E. Liatard, et al., Europhys. Lett. 13, 401 (1990).
  * (33) L. Nottale, Found Sci. 15, 101 (2010).


ssDNA short = ssDNA , long = single-stranded DNA dsDNA short = dsDNA , long =
double-stranded DNA GUV short = GUV , long = giant unilamellar vesicle TBE
short = TBE , long = Tris-borate-EDTA TE short = TE , long = Tris-EDTA DOPC
short = DOPC , long = 1,2-Dioleoyl-sn-glycero-3-phophocholine NBD-PE short =
NBD-PE , long = 1,2-dioleoyl-sn-
glycero-3-phosphoethanolamine-N-(7-nitro-2-1,3-benzoxadiazol-4-yl) (ammonium
salt) ITO short = ITO , long = indium tin oxide PDMS short = PDMS , long =
polydimethylsiloxane HPLC short = HPLC , long = high performance liquid
chromatography NaCl short = NaCl , long = sodium chloride

# Chemotactic crawling of multivalent vesicles along ligand-density gradients

Hannah Sleath  Department of Chemistry, Imperial College London, Molecular
Sciences Research Hub, 82 Wood Lane, London W12 0BZ, United Kingdom
fabriCELL, Imperial College London, Molecular Sciences Research Hub, 82 Wood
Lane, London W12 0BZ, United Kingdom  Bortolo Mognetti
<EMAIL_ADDRESS>Interdisciplinary Center for Nonlinear
Phenomena and Complex Systems, Université Libre de Bruxelles (ULB), B-1050
Brussels, Belgium  Yuval Elani<EMAIL_ADDRESS>Department of Chemical
Engineering, Imperial College London, Imperial College Road, London SW7 2AZ,
United Kingdom  fabriCELL, Imperial College London, Molecular Sciences
Research Hub, 82 Wood Lane, London W12 0BZ, United Kingdom  Lorenzo Di Michele
<EMAIL_ADDRESS>Department of Chemical Engineering and Biotechnology,
University of Cambridge, Philippa Fawcett Drive, Cambridge CB3 0AS, United
Kingdom  Department of Chemistry, Imperial College London, Molecular Sciences
Research Hub, 82 Wood Lane, London W12 0BZ, United Kingdom  fabriCELL,
Imperial College London, Molecular Sciences Research Hub, 82 Wood Lane, London
W12 0BZ, United Kingdom

###### Abstract

Living cells are capable of interacting with their environments in a variety
of ways, including cell signalling, adhesion, and directed motion. These
behaviours are often mediated by receptor molecules embedded in the cell
membrane, which bind specific ligands. Adhesion mediated by a large number of
weakly binding moieties—multivalent binding—is prevalent in a range of active
cellular processes, such as cell crawling and pathogen-host invasion. In these
circumstances, motion is often caused by gradients in ligand density, which
constitutes a simple example of chemotaxis. To unravel the biophysics of
chemotactic multivalent adhesion, we have designed an experimental system in
which artificial cell models based on lipid vesicles adhere to a substrate
through multivalent interactions, and perform chemotactic motion towards
higher ligand concentrations. Adhesion occurs _via_ vesicle-anchored receptors
and substrate-anchored ligands, both consisting of synthetic DNA linkers that
allow precise control over binding strength. Experimental data, rationalised
through numerical and theoretical models, reveal that motion directionality is
correlated to both binding strength and vesicle size. Besides providing
insights into the biophysics of chemotactic multivalent adhesion, our results
highlight design rules applicable to the development of biomimetic motile
systems for synthetic biology and therapeutic applications.

††preprint: APS/123-QED

##  I Introduction

Living cells interact with their environments _via_ receptors embedded in
their plasma membranes. These receptors can bind to specific ligand molecules,
resulting in processes vital to biological function such as adhesion to
external surfaces [1, 2] and the formation of pseudopodia involved in cell
motility [3]. Cell adhesion is often mediated by a large number of molecular
bonds between cell-membrane receptors and surface ligands. These multivalent
interactions can produce complex and useful emergent behaviours, such as
binding super-selectivity [4, 5], due to the interplay between enthalpic and
configurational effects [6]. Numerous experimental and theoretical studies
have been conducted to explore various aspects of multivalent interactions
between particles, membranes and surfaces, including the strength and rate of
adhesion [7, 8, 9, 10, 11]; the growth, size and stability of the contact
region between adhering objects [12, 13]; the self-assembly or fusion of
colloidal particles [14, 15, 16, 17, 18, 19, 20]; and receptor-mediated
endocytosis [21].

Cells and viruses adhering through multivalent interactions are known to
perform directed motion. Examples include the rolling and subsequent adhesion
of leukocytes to inflamed endothelia [22, 23]; and the diffusion of the
Influenza A [24, 25] and Herpes [26] viruses along the cell membrane prior to
invasion.

In many instances, motion directionality in multivalent adhesion systems is a
result of gradients in the type and density of ligands [27, 28, 29], making it
a form of chemotaxis conceptually analogous to the directional swimming of
algae, bacteria and sperm cells [30, 31, 32]. Despite relevant studies on
particle drift induced by shear flow [33, 34, 35, 36] and motion occurring in
low valency systems [37, 38], the biophysics of directed crawling mediated by
multivalent interactions remains largely unexplored. In particular, the
relative importance of factors such as binding strength and the size of the
adhering cell or particle on the ability to perform chemotactic crawling is
yet to be clarified.

Figure 1: Overview of the experimental system. (a) Schematic of the
experimental system, in which a giant unilamellar vesicle (GUV) functionalised
with DNA “receptor” constructs interacts with a surface density gradient of
complementary DNA “ligand” constructs attached to a substrate. Double-
cholesterol anchors are used for anchoring the constructs to the GUVs, while
biotin-streptavidin connections are used for the substrate. (b) Example of a
ligand-density gradient profile from experiments. (c) Microscopy images
showing the displacement of a GUV over 14 hours; the original position of the
vesicle is marked with a dashed circle. Axes indicating the directions
parallel and perpendicular to the gradient are included. (d) Mean displacement
over time of vesicles in systems with sticky end length $l=$ 5 nt. Data for
the mean displacement parallel and perpendicular to the ligand-density
gradient are shown with the standard error of the mean shaded. Positive values
of the displacement parallel to the gradient indicate motion towards higher
ligand-density regions. (e) Schematic of the theoretical model, illustrating
the motion of a GUV over a surface functionalised with ligands. The GUV is
attached to the surface by receptor-ligand bridges, which constrain its
position. Upon spontaneous unbinding of a bridge near its perimeter, the GUV
can explore beyond its original constrained region.

Here, we address these questions using a combination of theory, numerical
modelling, and experiments performed on a biomimetic model system. For
experiments, we consider synthetic cellular mimics consisting of GUV,
interacting with a solid surface through multivalent adhesion. As schematised
in Fig. 1a, both the “receptors” on the GUV and the “ligands” on the surface
consist of synthetic DNA constructs interacting through selective base-pairing
interactions that allow us to precisely modulate ligand-receptor affinity [7,
39, 13]. A gradient in the surface density of the ligands is established,
which causes the artificial cells to drift towards higher ligand
concentrations thanks to the reversibility of the ligand-receptor
interactions. We find that the vesicle drifting velocity is approximately
proportional to the unbinding rate of the ligand-receptor bridges, in
agreement with theoretical considerations. We further explore the relationship
between GUV size and motion, observing a positive correlation between vesicle
size and drifting velocity. Coarse-grained simulations based on the model
developed in Ref. [40] produce consistent trends in vesicle motion when
varying ligand-receptor affinity and vesicle size.

These results provide a deeper insight into the biophysics of multivalent
chemotactic motion, which could help to rationalise biological processes with
relevance to immunity, host-pathogen interactions and tissue dynamics.[22, 24,
27] Our work also provides design principles for engineering artificial cell
chemotaxis in the context of bottom-up synthetic biology, to generate
phenomena such as coordinated motion in multicellular systems, or targeting
and chasing of specific objects or chemical signals [41, 42, 43]. This could
lead to development of synthetic cellular solutions valuable in a variety of
applications _e.g._ therapeutics and targeted drug delivery [44, 45, 46, 47,
48, 49].

##  II Results and Discussion

###  II.1 Experimental system

A schematic of the experimental system is shown in Fig. 1a. Electroformed GUV
prepared from 1,2-Dioleoyl-sn-glycero-3-phosphocholine (DOPC), with typical
diameters of 1–50 $\mu$m (see Fig S1), are functionalised with DNA constructs,
here referred to as “receptors”. The receptors feature a double cholesterol
anchor that irreversibly partitions within the lipid bilayer. [7, 50] As DOPC
bilayers are fluid, receptors can freely diffuse laterally and redistribute
across the surface of the GUV. The surface density of receptors is
approximately 0.008 nm-2 (see Materials and Methods, Section IV). The
receptors interact with a second set of DNA constructs, here referred to as
“ligands”, which feature a biotin moiety for anchoring to a streptavidin-
coated substrate. The receptor and ligand constructs feature complementary
ssDNA sticky ends with lengths $l$ ranging from five to ten nucleotides (nt),
which can reversibly bind to each other. The constructs also feature rigid
dsDNA spacers designed to control their spatial extent, and short poly-T
domains to increase the flexibility and configurational freedom [7]. The DNA
sequences are listed in Tab. S1 with schematics of the constructs in Fig. S2.
We verified the correct assembly of the constructs and the binding of
complementary receptor-ligand pairs using agarose gel electrophoresis (see
Figs S3 and S4), while fluorescence microscopy was used to confirm attachment
of the receptor constructs to the vesicle membranes (see Fig. S5).

We set up a surface density gradient of ligands on the substrate, with density
typically varying from 0.01 nm-2 to 0.03 nm-2 over a distance of 100 $\mu$m
(see Fig. 1b). The method for generating the gradient is outlined in the
Materials and Methods (Section IV), and images of typical ligand-density
gradients are shown in Fig. S6. Receptor-functionalised GUV are then deposited
onto the gradient region where they adhere to the surface as sketched in Fig.
1a. For sufficiently strong adhesion, the GUV take the shape of a truncated
sphere, forming a flat adhesion patch visible in confocal cross sections (see
Fig. S7) [7, 39, 11]. Time lapse epifluorescence microscopy videos are
recorded, where both the GUV and the ligand-gradient can be visualised thanks
to calcein dye loaded in the vesicles and Alexa Fluor 488 modifications on a
sub-set of ligands. Example images in Fig. 1c show a GUV migrating in the
direction of the ligand-density gradient. Images are computationally segmented
to identify vesicle trajectories (see SI Section I). Projecting GUV
displacement onto the directions parallel and perpendicular to the local
gradient direction (see SI Section I C) allows us to quantify motion
directionality, as exemplified in Fig. 1d.

###  II.2 Numerical modelling

Coarse-grained simulations based on Ref. [40] are used to computationally
characterise the multivalent chemotactic system. We map the vesicles onto 2D
rigid disks of radius $R$, representing the perimeter of the flat adhesion
patch, as shown in Fig. 1e. In the following, we refer to the simulated
objects as vesicles or disks, interchangeably. Thermal fluctuations in vesicle
shape are neglected, and the disk-surface distance maintained constant. To
keep simulations affordable, we use disks with diameter $2R=0.2\,\mu$m,
1$\,\mu$m, and $1.8\,\mu$m, smaller by a factor $\sim 10$ compared to the GUV
used in experiments. The surface is randomly decorated by ligands to generate
a linear density profile with slope $\lambda$. The average surface density of
the ligands at the starting location of the disks is set to
$\rho_{\mathrm{L}}=0.021\,$nm-2, while the average density of receptors on the
vesicles is $\rho_{\mathrm{R}}=0.008\,$nm-2, consistently with the nominal
experimental values. The simulated gradient slope is set to
$\lambda=\rho_{\mathrm{L}}/10\,\mu$m-1, an order of magnitude higher than
experimental gradients such that the steepness of the gradient relative to
vesicle size is maintained. Note that further decreasing $\lambda$ would
rapidly make the system computationally untreatable, as for milder gradients
the drifting component of the motion is overwhelmed by the stochastic one.
Receptors on the vesicle can form bridges by binding free ligands within the
projection of the disk onto the surface. The rate constants at which bridges
form and break are given, respectively, by [7]

$\displaystyle
k_{\mathrm{on}}=\frac{\alpha\,k_{\mathrm{on}}^{\mathrm{sol}}}{3\pi R^{2}L}$
$\displaystyle
k_{\mathrm{off}}=\alpha\,\rho_{0}\,k_{\mathrm{on}}^{\mathrm{sol}}\exp[\beta\Delta
G_{0}],$ (1)

where $k_{\mathrm{on}}^{\mathrm{sol}}$ is the hybridisation rate of free
oligonucleotides in solution, $L$ is the length of the dsDNA spacer of the
ligands and receptors, $\Delta G_{0}$ is the standard hybridisation free
energy of the sticky ends, and $\rho_{0}=1$ M is the standard concentration.
$\alpha$ is a non-dimensional factor $<1$, accounting for the fact that the
hybridization kinetics are expected to be slower as a result of the DNA
constructs being tethered to surfaces. We set $\alpha=0.1$ as an estimate
based loosely on observations made in Ref. [5]. As simulated and theoretical
drift velocities are proportional to $\alpha$, changing this parameter affects
the absolute values but not the trends generated by varying other system
parameters. Vesicles are represented as hemispheres with a total surface area
equal to $3\pi R^{2}$. Accordingly, each free sticky end is taken as uniformly
distributed within a layer of thickness equal to the ligand/receptor length
($L$) surrounding the vesicle’s surface, and a resulting density equal to
$1/(3\pi R^{2}L)$. This density is used to calculate the binding rate
$k_{\mathrm{on}}$ using standard reaction equations [7]. By setting
$k_{\mathrm{on}}^{\mathrm{sol}}=10^{6}$ M-1s-1 (an approximate estimate for
short DNA oligomers [51]), $L=10$ nm, and $\Delta G_{0}$ values estimated
using the nearest neighbour thermodynamic model [52, 53, 54, 55], we obtain
$k_{\mathrm{off}}=$ 1.10 s-1, 0.174 s-1, and 0.0084 s-1 for sticky ends with
$l=$ 5, 6, and 7 nt, respectively.

In the model, we do not track the specific position of the receptors when free
(_i.e._ not bound to a ligand), consistently with the observation that the
receptors rearrange onto the bilayer much faster compared to the timescales of
vesicle motion. At any point in time, the set of configurations available to
the vesicle is limited by the formed bridges, as they are constrained to
remaining inside the perimeter of the disk. Until the bonds rearrange, the
disk can thus only rattle around the small sub-set of $n_{\mathrm{cb}}$
“constraining bridges” that can make direct contact with its perimeter (where
$\langle n_{\mathrm{cb}}\rangle=5$ [40]). Motion of the vesicle over larger
distances can only emerge following reconfiguration of the ligand-receptor
bridges, which thus limits motility [40].

We simulate the system using a reaction-diffusion algorithm [7, 56]. At each
simulation step, of duration $\Delta t$, we iteratively employ the Gillespie
algorithm [57] to break existing ligand-receptor bridges or form new ones, and
compute the time required for each reaction, until the total reaction time
exceeds $\Delta t$. Following the reaction step, we update the position of the
disk by randomly (uniformly) selecting a new centre of mass from all the
possible locations compatible with the new configuration of bridges [40]. The
simulation code can be found online [58].

###  II.3 Weaker binding enhances directional motion

Figure 2: Drift velocity increases with decreasing sticky end length. Data
illustrating the effect of sticky end length $l$ on the motion of the GUVs
along the ligand surface density gradient, for systems with $l=$ 5–7 nt.
Experimental data for vesicles of all sizes are shown in (a), while simulation
data for 1 $\mu$m-diameter vesicles are shown in (b). Solid lines indicate the
mean displacement of the vesicles parallel and perpendicular to the gradient
over time, where positive values of velocity parallel to the gradient indicate
motion towards regions of higher ligand density. Dashed lines indicated linear
fits to the data for parallel motion, with the slope annotated. Simulated
trajectories have been cropped to the same duration to enable visual
comparison, while the straight lines have been fitted to the entire, uncropped
trajectories (displayed in Figs S9 and S13). The scatter plots with marginal
histograms show the average velocities of individual vesicles parallel and
perpendicular to the gradient, calculated as total trajectory displacement
divided by duration. The y-axes of the histograms indicate the number of
vesicles in each bin. Red dashed lines and the red cross mark the mean of the
distributions. Note that the axes have been cropped to aid visualisation of
the majority of data points. Graphs including all data points are shown in
Figs S10 and S14.

A key factor affecting the strength of multivalent binding and the
corresponding vesicle mobility is the strength and (un)binding rates of
individual receptor-ligand interactions. This factor can be tuned by varying
the length of the complementary sticky ends; in our experiments and
simulations we studied systems with sticky end lengths $l=$ 5, 6 and 7 nt.

Experimentally, under all these conditions, the vesicles adhere to the
substrate and form a stable, flat adhesion patch, as confirmed by confocal
cross sections (see Fig. S7). As discussed in Section I D of the SI, in all
samples, we notice that a small proportion of vesicles (approximately 5–6%)
remains much more static that the rest of the population, indicated by a
significantly smaller diffusion coefficient (see Fig. S8). We believe that
this effectively immobile sub-population could emerge due to a number of
factors, such as sample impurities or surface defects being erroneously
identified as vesicles; non-specific adhesion or trapping of the vesicles on
the substrate; or tracking inaccuracies. In the following data analysis the
immobile population is excluded, although all plots in the main text have been
replicated in the SI with the immobile population included for transparency
(Figs S9–14), showing no qualitative changes to the findings outlined below.

To study the speed and direction of the vesicle motion in both experiments and
simulations, each trajectory was separated into components parallel and
perpendicular to the local density gradient of DNA ligand constructs.
Experimental and simulated data are collated in Fig. 2a and 2b. Here we report
the mean displacement projected along the directions parallel and
perpendicular to the ligand-density gradient (left) and the two-dimensional
distributions of average velocities, as computed from the initial and final
positions of the vesicles in the trajectories (right). From the average
velocity distributions, we note that, in all cases, the motion is
predominantly stochastic. However, while the projection of the average
velocity onto the direction perpendicular to the gradient is centered around
zero, _i.e._ , shows no directional bias, the projection along the direction
of the gradient has non-zero mean for $l=$ 5 and 6 nt, both in experiments and
simulations. The observed bias indicates motion towards denser regions of the
ligand carpet, as expected. For $l=7$ nt, no clear directional bias is noted,
in simulations or experiments. The width of the average velocity distributions
also decreases with increasing $l$, indicating that vesicles with shorter
sticky ends diffuse less, regardless of directionality, as shown in Fig. S15.
This relative decrease is more pronounced in simulations than it is in
experiments, as expected given that experimental data are inevitably impacted
by static localisation errors.

Directional bias of the motion towards higher ligand densities is better
visualised from the mean displacement data shown as a function of time (Figs
2a and 2b, left). In all cases, data on the motion perpendicular to the
gradient average to zero (within standard error), while data on the parallel
motion display a non-zero slope for $l=5$ and 6 nt, with the latter being less
pronounced. For $l=7$ nt, mean displacement parallel to the gradient shows a
negligible slope.

From the scatter plots in Fig. 2 we calculate the average experimental drift
velocities for the $l=$ 5 and 6 nt systems to be $5.3\times 10^{-3}$ nm s-1
and $1.2\times 10^{-3}$ nm s-1 respectively. Over experimental timescales
(_i.e._ 10 hours), this equates to total displacements of 190 nm and 43 nm
respectively, which are several orders of magnitude smaller than the typical
vesicle diameter ($\approx 10$ $\mu$m; see Fig. S1). However, we note that the
distribution of vesicle displacements parallel to the gradient has a positive
skew for both $l=$ 5 and 6 nt, with maximum displacements along the gradient
of 5.8 $\mu$m and 5.1 $\mu$m respectively, comparable with the average vesicle
radius. For comparison, a freely diffusing vesicle of this size would travel
$\approx$ 80 $\mu$m along the surface in this time scale (using the Stokes-
Einstein equation to estimate the diffusion constant), which shows that the
vesicle mobility is reduced by the multivalent adhesion.

To put the extent of the vesicle displacements into perspective, we consider
the spacing between ligand-receptor bridges. In Fig. S16 we see that the
system is in a saturation regime where the ligand-receptor bridge density is
limited by either the density of ligands or the density of receptors. This
means that the bridge density is on the order of 0.01 nm-2, which corresponds
to an inter-bridge spacing of 10 nm. A vesicle displacement of $\sim 5\,\mu$m
is several orders of magnitude higher than the inter-bridge spacing,
indicating that the observed motion cannot simply be due to the rattling of
vesicles around fixed bridging points but must require the formation and
breakage of a large number of bridges.

Furthermore, we note that a typical vesicle with a diameter of 10 $\mu$m,
located centrally in the gradient region of Fig. 1b, spans ligand densities
that only differ by $\sim 10$% (from 0.019 nm-2 to 0.021 nm-2). Therefore,
during each formation of a new ligand-receptor bridge, it is only 10% more
likely (at most) for the bridge to form on the side of the contact region
facing towards higher ligand density, compared to the opposite side. This
small difference in binding probability across the span of the vesicle
suggests that a large number of bridge formation and breakage events are
required for the directional vesicle drift along the gradient to be
distinguishable from stochastic motion.

When comparing experiments with simulated trajectories, we note that the
average drift velocities from experimental data for $l=$ 5 and 6 nt are an
order of magnitude higher than the simulated velocities (see Fig. 2). This
deviation could be due to the simulated vesicle diameters being an order of
magnitude smaller than in experiments, although the simulated gradients have
been adjusted to achieve the same gradient steepness relative to vesicle size.
The effect of vesicle size on drift velocity is further discussed in Section
II.4. The difference in drift velocities between experiments and simulations
could also be due to differences in the gradient profiles or to small
convective flows and tracking inaccuracies inflating experimental estimates.
Furthermore, calculations of the bridge formation and breakage rates (Eq. 1)
for the simulations involve assumption of the correction factor $\alpha$,
which could introduce further discrepancies.

Besides systems with sticky end lengths $l=$ 5–7 nt, we also experimentally
tested samples with shorter sticky end lengths ($l=$ 3 and 4 nt), with the
goal of further increasing vesicle mobility. However, we found that the
binding strength was insufficient to keep the GUV adhered to the surface,
causing them to drift driven by convection or gravity (see SI Fig. S17). We
additionally designed constructs with longer sticky ends ($l=$ 8–10 nt), which
showed negligible drift consistent with the data in Fig. 2 for $l=7$ nt.

####  II.3.1 Comparison of experimental drift velocities with theory

For a given number of ligands and receptors, the drift velocity is expected to
be a function of the off rate ($k_{\mathrm{off}}$), the number of bridges
($n_{\mathrm{b}}$), and the diffusion constant of free (non-adhering) vesicles
($D$). In systems featuring many bridges or low unbinding rates, the
displacement dynamics of the vesicles is limited by the timescales of
ligand/receptor reactions. This can be verified by comparing two
characteristic timescales. The first is the typical timescale through which
the translational configurational space available to the vesicles changes
following a binding and unbinding of ligands/receptors, which can be estimated
as $\tau_{\mathrm{r}}=1/(2k_{\mathrm{off}}n_{\mathrm{cb}})$, with $\langle
n_{\mathrm{cb}}\rangle\approx 5$ [40]. The second is the diffusion timescale
needed by a vesicle to explore the available space: $\tau_{\mathrm{d}}=\pi
R^{2}n_{\mathrm{cb}}/(Dn_{\mathrm{b}}^{2})$.

In conditions considered here, $\tau_{\mathrm{r}}$ is at least 7 orders of
magnitude bigger than $\tau_{\mathrm{d}}$ (see SI Section II A), making the
dynamics of our systems reaction-limited. Under these circumstances, we expect
$v_{\mathrm{drift}}$ to be entirely determined by $k_{\mathrm{off}}$,
resulting in the following scaling

$v_{\mathrm{drift}}\propto k_{\mathrm{off}}f(R,n_{\mathrm{b}}),$ (2)

where $f(R,n_{\mathrm{b}})$ is a general function on the variables $R$ and
$n_{\mathrm{b}}$.

For comparison, we derive an analytical prediction of the drift velocity. In
our previous contribution, we predicted that a receptor-decorated disk on a
uniformly ligand-decorated surface would diffuse with $D\sim
k_{\mathrm{off}}R^{2}/\langle n_{\mathrm{b}}\rangle^{2}$ [40]. By applying the
fluctuation-dissipation theorem to this expression, we find (see SI Section II
B)

$v_{\mathrm{drift}}^{\mathrm{FD}}=-\beta
D\frac{\mathrm{d}V_{\mathrm{multi}}}{\mathrm{d}x}=-\frac{\lambda
k_{\mathrm{off}}}{2\langle\rho_{\mathrm{b}}\rangle^{2}}\log\left[1-\frac{\langle\rho_{\mathrm{b}}\rangle}{\langle\rho\rangle}\right],$
(3)

where $V_{\mathrm{multi}}$ is the multivalent free-energy [59, 7] and
$\langle\rho\rangle$ ($\langle\rho_{\mathrm{b}}\rangle$) the density of
(bound) ligands found at the center of the disk. The expression in Eq. 3
agrees with the relationship anticipated in Eq. 2. However, the theoretical
prediction in Eq. 3, substantially overestimate drift velocities determined in
both experiments and simulations, as shown in Figs S18 and S19. The
discrepancies between Eq. 3 and simulations decrease at low values of
$n_{\mathrm{b}}$. This discrepancy can be rationalised by noting that, for
large $n_{\mathrm{b}}$ (and thus large drift forces), our systems may deviate
from the linear response regime under which the fluctuation-dissipation
theorem applies. This observation is in agreement with previous findings that
coarse-grained dynamics in non-equilibrium conditions are not necessarily
described by a simple Langevin equation [60, 61].

Although the theoretical prediction in Eq. 3 does not appear to be valid under
experimentally relevant conditions, we still expect the drift velocity to
follow the general scaling of Eq. 2. As previously mentioned, Fig. S16
indicates that for $l=$ 5,6 and 7 nt, the system is in a saturation regime
where the bridge density $\langle n_{\mathrm{b}}\rangle$ is limited by either
the ligand density or the receptor density. It follows that the $R$ and
$n_{\mathrm{b}}$ dependencies in Eq. 2 should be unchanged when comparing
experiments with different $l$, resulting in $v_{\mathrm{drift}}\propto
k_{\mathrm{off}}$.

To verify the predicted proportionality, we compute the ratio between
$v_{\mathrm{drift}}$ values for systems with different $l$, and compare it
with the ratio between $k_{\mathrm{off}}$ values. From the linear fits in Fig.
2, we extract $v^{l=5}_{\mathrm{drift}}/v^{l=6}_{\mathrm{drift}}=5.84$ and
4.90 for experiments and simulations, respectively. Both these values are in
good agreement with $k^{l=5}_{\mathrm{off}}/k^{l=6}_{\mathrm{off}}=6.30$,
obtained using the values reported in Section II.2. The small discrepancy
possibly derives from using the same $\alpha$ for both sticky ends when
computing $k_{\mathrm{off}}$ (Eq. 1), with $\alpha$ being, in principle,
sequence dependent. While similar comparisons could not be carried out for the
low-mobility system with $l=7$ nt, the analysis for $l=5$ nt and $l=6$ nt
supports a scenario where the drift velocity is dictated by unbinding rates,
in alignment with the model leading to Eq. 2.

###  II.4 Larger GUV exhibit chemotaxis to a greater extent

We also experimentally observe that the size of the GUV has an effect on the
drift velocity along the gradient. When binning the experimental trajectories
by vesicle diameter, we notice that larger vesicles travel on average a
greater distance along the gradient, as shown in Fig. 3 (left). The trend is
most evident in the $l=5$ nt system, where vesicles are the most mobile, but
also present for the less mobile $l=6$ nt system. Simulated trajectories for
vesicles of varying diameters confirm the experimental trends for $l=5$ and 6
nt (see Fig. 3, right). In both the experiments and simulations, we found that
the average vesicle drift velocity is typically greater than the median, due
to the velocity distribution having a positive skew. We do not observe a clear
relationship between vesicle size and motion along the gradient for $l=7$ nt,
in neither experiments nor simulations, likely due to any signal being drowned
by noise in this low-mobility system. Regardless of size, we observe that, on
average, vesicles drift a negligible distance along the direction
perpendicular to the gradient, as shown in Fig. S12 and S21 and consistent
with the data in Fig. 2.

Figure 3: Drift velocity increases with increasing vesicle diameter. Box plots
illustrating the relationship between GUV size and the average drift velocity
of along the ligand density gradient, for sticky length $l=$ 5, 6 and 7 nt.
Data show the distribution of average velocity along the gradient, which has
been calculated as total displacement parallel to the gradient divided by
trajectory duration. Positive values indicate motion of the towards regions of
higher ligand surface density. For each box plot, the mean is marked with a
cross. Note that outliers have been excluded for better visualisation. Graphs
with outliers included are presented in Fig. S10 and S14.

One possible explanation for the relationship between vesicle size and drift
velocity in the mobile systems ($l=$ 5 and 6 nt) is that larger vesicles cover
a greater extent of the ligand density gradient, corresponding to a greater
change in binding affinity across the width of the vesicle. Following from the
discussion in Section II.3, for a vesicle located centrally on the gradient
depicted in Fig. 1b, the ligand density varies by $\sim$5% across a 5
$\mu$m-diameter contact region, and by $\sim$20% across a 20 $\mu$m-diameter
contact region. This means that the larger vesicle would have a higher rate of
bridge formation on the side of the contact region facing up the gradient
compared with the smaller vesicle, resulting in a greater drift velocity. This
effect has also been discussed in Ref. [36], where the motility of influenza
virus particles adhered to surface-bound molecular density gradients was
studied. In this case, no directional bias was observed in the virus motion,
which the authors ascribed to the small size of the virus particles relative
to the steepness of the ligand-density gradient.

Finally, we note that while larger vesicles can more effectively sense
gradients, they are also less mobile due to the larger number of bridges [40].
These two competing factors perfectly compensate each other in the theoretical
prediction reported in Eq. 3, where $v^{\mathrm{FD}}_{\mathrm{drift}}$ is not
dependent on $R$. Evidence of a size-dependent drift velocity in both
experiments and simulations further highlights the limitations of the
theoretical approach leading to Eq. 3 in the strong-binding regime relevant to
our system.

##  III Conclusions

In summary, through experiments, theory and simulations, we studied the
directional motion of receptor-decorated lipid vesicles adhering to a ligand-
decorated surface, in the presence of a gradient in ligand density – a model
system for chemotactic crawling. Experimentally, both receptors and ligands
consisted of synthetic DNA constructs, the former connected to giant
unilamellar lipid vesicles and the latter to a solid surface. Coarse-grained
simulations relied on a multi-scale approach first reported in Ref. [40],
which accurately describes the binding and unbinding dynamics of ligands and
receptors through a Gillespie algorithm.

Both simulations and experiments showed directional motion of the vesicles
towards ligand-dense regions of the surface. The magnitude of the directional
drift was observed to decrease with increasing ligand-receptor affinity, which
could be easily controlled by changing the length of the single-stranded DNA
sticky ends through which our constructs interact.

With theoretical arguments we demonstrated that, in the regime relevant to our
system, directional motion is limited by the timescales of the ligand-receptor
reactions. The corresponding relationship between drift velocity and ligand-
receptor unbinding rates were confirmed by both experiments and simulations.

We also observed a positive correlation between vesicle size and drift
velocity in both experiments and simulations. This trend, more noticeable in
more mobile systems, is explained by the ability of larger vesicles to probe a
greater extent of the ligand density gradient, thus generating greater spatial
asymmetry in the ligand-receptor binding probabilities.

Our findings offer quantitative insights on the mechanisms underpinning
directional motion in multivalent systems, highly relevant to a variety of
biological processes involving the membrane interactions. Examples include
immune cells adhesion [22, 23], viral invasion [24, 25, 26], and tissue
dynamics [22, 24, 27]. In future studies, it would be interesting to explore
parameter spaces such as temperature, buffer conditions and gradient
steepness, as well as working in more weakly-binding regimes (_i.e._ by
reducing sticky end length or decreasing DNA density); this may require
alteration of the experimental setup to reduce convection effects.

Our models and experimental implementation will also support efforts to
engineer directed motion in synthetic cellular systems. Mobile synthetic cells
capable of performing simple chemotaxis would indeed be valuable for a vast
range of applications, such as: smart drug delivery systems that can target
and release drugs at specific locations in the body [44, 45, 46, 47, 48, 49];
vesicle-based biosensors or biomedical imaging systems that signal the
presence of specific molecules or conditions [62, 63, 64]; or tissue
engineering applications where gradients of molecules or growth factors are
used to direct cell migration and tissue formation [65, 66, 67].

##  IV Materials and Methods

###  IV.1 Experimental materials

##### DNA oligonucleotides.

Sequences of the DNA oligonucleotides were designed using the NUPACK design
tool [68]. DNA strands with 5’-cholesteryl modifications were purchased from
Eurogentec (Liege, Belgium) with HPLC purification. All other DNA strands were
purchased from Integrated DNA Technologies (Coralville, Iowa, United States).
Modified strands were purified via HPLC, and unmodified strands were purified
by standard de-salting.

##### Buffers.

NaCl (BioUltra, $>$99.5$\%$) and 100$\times$ TE buffer were purchased from
Sigma-Aldrich (Gillingham, UK). 10$\times$ TBE was purchased from ThermoFisher
Scientific (UK). All buffer solutions were diluted with the appropriate amount
of MilliQ water, and were filtered through 0.22 $\mu$m pore poly(ether
sulfone) filters (Millex) prior to use.

##### Gel electrophoresis.

Agarose was purchased from Sigma-Aldrich (Gillingham, UK). Ultra Low Range DNA
Ladder and 6$\times$ TrackIt Cyan/Yellow Loading Buffer were purchased from
ThermoFisher Scientific (UK). 10,000$\times$ SYBR Safe DNA gel stain was
purchased from APExBIO Technology LLC (Houston, Texas, USA).

##### Vesicle generation.

D-(+)-glucose ($>$99.5$\%$), sucrose, chloroform and indium tin oxide coated
glass slides (surface resistivity 15-25 $\Omega$ sq-1) were purchased from
Sigma-Aldrich (Gillingham, UK). $>$99$\%$ DOPC and NBD-PE was purchased from
Avanti Polar Lipids Inc. (Alabaster, Alabama, USA).

##### Well plates and adhesive covers.

Streptavidin-coated high capacity 96-well strip plates were purchased from
Sigma-Aldrich (Gillingham, UK). Adhesive clear foils for 96-well plates were
purchased from Sarstedt Ltd (Leicestershire, UK).

###  IV.2 Experimental methods

##### DNA Reconstitution.

The DNA strands were shipped lyophilised and were reconstituted to a
concentration of approximately 100 $\mu$M in 1$\times$ TE buffer (1 mM EDTA,
10 mM Tris, pH 8.0). The concentration of the reconstituted DNA was determined
via UV-vis spectrophotometry on a Thermo Scientific NanoDrop One, by measuring
the ratio of absorbance at 260 nm to the sequence-specific extinction
coefficient of the DNA. Stock solutions of DNA were then stored at
$-$20$\degree$C.

##### DNA Hybridisation.

Assembly of the multistranded receptor and ligand constructs was facilitated
by thermal annealing. For each construct, the constituent strands were diluted
to 2 $\mu$M and combined in TE buffer containing 100 mM NaCl. The strands were
first heated to 90$\degree$C for 5 minutes to ensure melting of all DNA
duplexes, and then cooled to 20$\degree$C at a rate of $-$0.5$\degree$C min-1,
on a Bio Rad C1000 Thermal Cycler. The sizes of the annealed DNA constructs
were evaluated _via_ agarose gel electrophoresis to verify correct folding, as
reported in Figs S3 and S4.

##### Agarose gel electrophoresis.

3% (w/v) agarose gels were casted in 1$\times$ TBE buffer with 1$\times$ SYBR
Safe DNA gel stain. The 10 $\mu$l wells were loaded with approximately 300 ng
DNA samples with loading dye; the outer two wells were loaded with DNA
ladders, to enable comparison with the DNA samples. The gels were run for 90
minutes at 120 V (electric field strength of 6 V cm-1) and then imaged using a
Syngene Dyversity 4 gel imager.

##### Microscopy.

Epifluorescence microscopy imaging was performed on a Nikon Eclipse Ti2-E
inverted microscope using a Nikon CFI Plan Apochromat Lambda D 10$\times$ dry
objective (NA 0.45). Confocal microscopy imaging was performed on a Leica TCS
SP5 Confocal microscope using a HCX PL APO CS 63.0$\times$ oil-immersion
objective (NA 1.40).

##### Preparing vesicles.

The GUV were produced _via_ electroformation, a straightforward and commonly
used method for vesicle generation [69]. Lipid solutions consisting of DOPC
with 1 mol$\%$ fluorescent lipid NBD-PE were prepared by dissolving the
appropriate quantities of lipids in chloroform to yield a 1 mg ml-1 solution.
30 $\mu$l of this solution was spread evenly on an ITO slide and vacuum-
desiccated for 30 minutes to remove residual chloroform, resulting in a lipid
film. A 5 mm-thick PDMS spacer with a central cut-out was sandwiched between
this film and another ITO slide to create a chamber, and the chamber was
filled with a solution consisting of 300 mM sucrose and 50 $\mu$M calcein in
water. A function generator (Aim-TTi, TG315, Huntingdon, UK) was used to apply
an alternating electric field at a peak-to-peak voltage of 1.5 V was applied
across the ITO slides at 10 Hz for two hours, followed by one further hour at
2 Hz. During this time, the chamber was left in a 60$\degree$C oven to ensure
that the lipid was in the fluid phase. Finally, the chamber was opened and the
resulting vesicles were collected. The vesicle sample was inspected _via_
epifluorescence microscopy, by exciting both the NBD-PE fluorophore and the
encapsulated calcein. Note that the calcein signal overwhelms the signal from
the NBD-PE; the latter was included to aid visualisation of the vesicles in
the case of any calcein leakage. Ease of imaging can be improved by mixing the
vesicle sample with a solution of 300 mM glucose in MilliQ water, such as in a
1:10 ratio of vesicles to glucose solution. This introduces a density
difference which causes the vesicles to settle at the bottom of the viewing
chamber, and a refractive index contrast to enable viewing in brightfield. The
size distributions of vesicles in our experiments are shown in Fig. S5.

##### Vesicle functionalisation.

GUV were functionalised by incubation with cholesterolised DNA constructs. A
sample was prepared with 20 v/v$\%$ 2 $\mu$M DNA constructs in 100 mM NaCl TE
buffer, 11 v/v$\%$ GUV (prepared as detailed in the previous paragraph) in 300
mM glucose solution, and 69 v/v$\%$ “correction buffer”. This correction
buffer consists of 116 mM NaCl, 78.2 mM glucose and 1.16$\times$ TE buffer,
and was included to match the osmolarity of the interior and exterior of the
GUV and thus to prevent osmotic-shock-induced rupture. The sample was left to
incubate on a roller mixer overnight. After the functionalisation process, the
vesicle sample was washed to minimise the amount of DNA constructs remaining
in solution that have not been incorporated into the lipid membranes. This was
achieved by leaving the sample to stand for at least 15 minutes to allow the
vesicles to settle, and then replacing approximately 90$\%$ of the eluent with
an iso-osmolar solution consisting of 100 mM NaCl and 87 mM glucose in TE
buffer. The washing process was carried out at least three times. To check
that the vesicles had been successfully functionalised and the lipid membranes
were saturated with DNA constructs, the vesicles were prepared with
fluorescently-tagged lipids and DNA constructs, and imaged _via_
epifluorescence and confocal microscopy. The estimated density of DNA
constructs in the lipid membrane is approximately 0.008 nm-2.

##### Ligand-density gradient generation and characterisation.

A 0.5 $\mu$l droplet of 2 $\mu$M biotinylated DNA “ligand” constructs in
buffer solution (100 mM NaCl and 1$\times$TE) was deposited on the substrate,
a well plate pre-coated with streptavidin. This was left to incubate for five
minutes to allow the DNA constructs to bind to the substrate. Excess DNA was
then washed away by quickly rinsing the substrate with buffer solution (100 mM
NaCl and 1$\times$TE). This left a coating of DNA in the circular region that
had been covered by the droplet, with a diameter of $\approx$3.5 mm and a
ligand density of $\approx$0.033 nm-2. Due to the fast rate of binding between
streptavidin and biotin, a small amount of excess DNA bound to the surrounding
region during the washing process; the surface density of DNA in this region
was measured to be 0.005–0.01 nm-2 on average. At the perimeter of the
circular maximally-coated region, there was a steep surface density gradient
spanning a distance of approximately 300 $\mu$m. The ligand constructs were
tagged with Alexa Fluor 488, enabling the ligand coverage to be imaged and
characterised by excitation of the fluorophore. The ligand density in the
maximum-density region of the wells was estimated by applying a droplet of
ligand solution to the well for five minutes, like in the gradient generation
procedure; the supernatant was diluted and its fluorescence was recorded on a
BMG LABTECH CLARIOstar Plus microplate reader. The fluorescence of the diluted
supernatant was compared with the fluorescence of ligand solutions of known
concentrations, to estimate the quantity of ligand constructs that had been
deposited on the substrate. To characterise the ligand-density gradient
profiles in experiments, the ligand coverage was imaged with epifluorescence
microscopy. Due to non-uniform illumination of the sample, the fluorescence
images were normalised by a control image taken within a non-functionalised
well containing calcein solution (see Fig. S6). The ligand density was then
estimated by interpolating between the background fluorescence signal
(recorded from a non-functionalised well) and the average fluorescence signal
of a maximum ligand-density region of the well.

##### Recording vesicle motion.

After generating the ligand-density gradients, the receptor-coated vesicles
were added to the wells, and the well plate was covered with a transparent
adhesive foil to minimise evaporation and evaporation-driven flows in the
sample. The sample was then left to sit for one hour to allow most of the
vesicles to sink to the bottom of the well and adhere. A time lapse recording
was then taken of the vesicles on the substrate, over a duration of 10–15
hours.

##### Vesicle tracking

The time lapse recordings were initially stabilised using the Image Stabilizer
plugin on ImageJ [70] to remove any unwanted effects of stage drift or
inaccuracies in stage positioning. The image data were then loaded into
Python, and processed to remove noise and prepare the images for the tracking
algorithm. Vesicles were located in each frame of the time lapse recordings
using a template-matching method. Particle trajectories were then identified
by linking particle positions between successive time frames; this was
achieved using the Python module Laptrack [71, 72]. More details on data
processing and the vesicle locating and tracking algorithms can be found in
Section I of the SI.

###### Acknowledgements.

H.S. acknowledges funding from an EPSRC DTP Studentship in the Department of
Chemistry at Imperial College London (EP/R513052/1). L.D.M. acknowledges
support from the European Research Council (ERC) under the Horizon 2020
Research and Innovation Programme (ERC-STG No 851667 – NANOCELL) and a Royal
Society University Research Fellowship (UF160152, URF\R\221009). B.M.M. is
supported by a PDR grant of the FRS-FNRS (Grant No. T.0158.21). Y.E. is
supported by a UK Research and Innovation Future Leaders Fellowship
(MR/S031537/1).

## References

  * Ahmad Khalili and Ahmad [2015] A. Ahmad Khalili and M. R. Ahmad, International journal of molecular sciences 16, 18149 (2015).
  * Kam and Boxer [2001] L. Kam and S. G. Boxer, Journal of Biomedical Materials Research: An Official Journal of The Society for Biomaterials, The Japanese Society for Biomaterials, and The Australian Society for Biomaterials and the Korean Society for Biomaterials 55, 487 (2001).
  * Van Haastert and Devreotes [2004] P. J. Van Haastert and P. N. Devreotes, Nature reviews Molecular cell biology 5, 626 (2004).
  * Liu  _et al._ [2020] M. Liu, A. Apriceno, M. Sipin, E. Scarpa, L. Rodriguez-Arco, A. Poma, G. Marchello, G. Battaglia, and S. Angioletti-Uberti, Nature communications 11, 4836 (2020).
  * Linne  _et al._ [2021] C. Linne, D. Visco, S. Angioletti-Uberti, L. Laan, and D. J. Kraft, Proceedings of the National Academy of Sciences 118, e2106036118 (2021).
  * Varilly  _et al._ [2012] P. Varilly, S. Angioletti-Uberti, B. M. Mognetti, and D. Frenkel, The Journal of chemical physics 137, 094108 (2012).
  * Mognetti  _et al._ [2019] B. M. Mognetti, P. Cicuta, and L. Di Michele, Reports on progress in physics 82, 116601 (2019).
  * Angioletti-Uberti [2017] S. Angioletti-Uberti, npj Computational Materials 3, 48 (2017).
  * Curk  _et al._ [2018] T. Curk, J. Dobnikar, and D. Frenkel, Multivalency: Concepts, Research & Applications , 75 (2018).
  * Lanfranco  _et al._ [2019] R. Lanfranco, P. K. Jana, L. Tunesi, P. Cicuta, B. M. Mognetti, L. Di Michele, and G. Bruylants, Langmuir 35, 2002 (2019).
  * Amjad  _et al._ [2017] O. A. Amjad, B. M. Mognetti, P. Cicuta, and L. Di Michele, Langmuir 33, 1139 (2017).
  * Shenoy and Freund [2005] V. Shenoy and L. Freund, Proceedings of the National Academy of Sciences 102, 3213 (2005).
  * Parolini  _et al._ [2015] L. Parolini, B. M. Mognetti, J. Kotar, E. Eiser, P. Cicuta, and L. Di Michele, Nature communications 6, 5948 (2015).
  * Beales and Vanderlick [2007] P. A. Beales and T. K. Vanderlick, The Journal of Physical Chemistry A 111, 12372 (2007).
  * Beales  _et al._ [2011] P. A. Beales, J. Nam, and T. K. Vanderlick, Soft Matter 7, 1747 (2011).
  * Beales and Vanderlick [2014] P. A. Beales and T. K. Vanderlick, Advances in colloid and interface science 207, 290 (2014).
  * Rogers  _et al._ [2016] W. B. Rogers, W. M. Shih, and V. N. Manoharan, Nature Reviews Materials 1, 1 (2016).
  * Nykypanchuk  _et al._ [2008] D. Nykypanchuk, M. M. Maye, D. Van Der Lelie, and O. Gang, Nature 451, 549 (2008).
  * Flavier and Boxer [2017] K. M. Flavier and S. G. Boxer, Biophysical Journal 113, 1260 (2017).
  * Chan  _et al._ [2009] Y.-H. M. Chan, B. van Lengerich, and S. G. Boxer, Proceedings of the National Academy of Sciences 106, 979 (2009).
  * Di Michele  _et al._ [2018] L. Di Michele, P. K. Jana, and B. M. Mognetti, Physical Review E 98, 032406 (2018).
  * Ley  _et al._ [2007] K. Ley, C. Laudanna, M. I. Cybulsky, and S. Nourshargh, Nature Reviews Immunology 7, 678 (2007).
  * Muller [2002] W. A. Muller, Laboratory investigation 82, 521 (2002).
  * de Vries  _et al._ [2020] E. de Vries, W. Du, H. Guo, and C. A. de Haan, Trends in microbiology 28, 57 (2020).
  * Vahey and Fletcher [2019] M. D. Vahey and D. A. Fletcher, Elife 8, e43764 (2019).
  * Delguste  _et al._ [2018] M. Delguste, C. Zeippen, B. Machiels, J. Mast, L. Gillet, and D. Alsteens, Science advances 4, eaat1273 (2018).
  * Basara  _et al._ [1985] M. L. Basara, J. B. McCarthy, D. W. Barnes, and L. T. Furcht, Cancer research 45, 2487 (1985).
  * Thibault  _et al._ [2007] M. M. Thibault, C. D. Hoemann, and M. D. Buschmann, Stem cells and development 16, 489 (2007).
  * Haessler  _et al._ [2011] U. Haessler, M. Pisano, M. Wu, and M. A. Swartz, Proceedings of the National Academy of Sciences 108, 5614 (2011).
  * Adler [1966] J. Adler, Science 153, 708 (1966).
  * Choi  _et al._ [2016] H. I. Choi, J. Y. H. Kim, H. S. Kwak, Y. J. Sung, and S. J. Sim, Biomicrofluidics 10, 014121 (2016).
  * Eisenbach and Giojalas [2006] M. Eisenbach and L. C. Giojalas, Nature reviews Molecular cell biology 7, 276 (2006).
  * Hammer and Lauffenburger [1987] D. A. Hammer and D. A. Lauffenburger, Biophysical journal 52, 475 (1987).
  * Dasanna  _et al._ [2017] A. K. Dasanna, C. Lansche, M. Lanzer, and U. S. Schwarz, Biophysical journal 112, 1908 (2017).
  * Porter  _et al._ [2021] C. L. Porter, S. L. Diamond, T. Sinno, and J. C. Crocker, Biophysical Journal 120, 2102 (2021).
  * Hamming  _et al._ [2023] P. E. Hamming, N. J. Overeem, K. Diestelhorst, T. Fiers, M. Tieke, G. M. Vos, G.-J. P. Boons, E. van der Vries, S. Block, and J. Huskens, ACS applied materials & interfaces (2023).
  * Perl  _et al._ [2011] A. Perl, A. Gomez-Casado, D. Thompson, H. H. Dam, P. Jonkheijm, D. N. Reinhoudt, and J. Huskens, Nature chemistry 3, 317 (2011).
  * Chang  _et al._ [2007] T. Chang, D. I. Rozkiewicz, B. J. Ravoo, E. Meijer, and D. N. Reinhoudt, Nano letters 7, 978 (2007).
  * Shimobayashi  _et al._ [2015] S. Shimobayashi, B. M. Mognetti, L. Parolini, D. Orsi, P. Cicuta, and L. Di Michele, Physical Chemistry Chemical Physics 17, 15615 (2015).
  * Lowensohn  _et al._ [2022] J. Lowensohn, L. Stevens, D. Goldstein, and B. M. Mognetti, The Journal of Chemical Physics 156, 164902 (2022).
  * Blain and Szostak [2014] J. C. Blain and J. W. Szostak, Annual review of biochemistry 83, 615 (2014).
  * Guindani  _et al._ [2022] C. Guindani, L. C. da Silva, S. Cao, T. Ivanov, and K. Landfester, Angewandte Chemie 134, e202110855 (2022).
  * Siton-Mendelson and Bernheim-Groswasser [2016] O. Siton-Mendelson and A. Bernheim-Groswasser, Cell adhesion & migration 10, 461 (2016).
  * Fischbach  _et al._ [2013] M. A. Fischbach, J. A. Bluestone, and W. A. Lim, Science translational medicine 5, 179ps7 (2013).
  * Joseph  _et al._ [2017] A. Joseph, C. Contini, D. Cecchin, S. Nyberg, L. Ruiz-Perez, J. Gaitzsch, G. Fullstone, X. Tian, J. Azizi, J. Preston, _et al._ , Science Advances 3, e1700362 (2017).
  * Kershaw  _et al._ [2002] M. H. Kershaw, G. Wang, J. A. Westwood, R. K. Pachynski, H. L. Tiffany, F. M. Marincola, E. Wang, H. A. Young, P. M. Murphy, and P. Hwu, Human gene therapy 13, 1971 (2002).
  * Moon  _et al._ [2011] E. K. Moon, C. Carpenito, J. Sun, L.-C. S. Wang, V. Kapoor, J. Predina, D. J. Powell Jr, J. L. Riley, C. H. June, and S. M. Albelda, Clinical cancer research 17, 4719 (2011).
  * Pockaj  _et al._ [1994] B. A. Pockaj, R. M. Sherry, J. P. Wei, J. R. Yannelli, C. S. Carter, S. F. Leitman, J. A. Carasquillo, S. M. Steinberg, S. A. Rosenberg, and J. C. Yang, Cancer 73, 1731 (1994).
  * Karp and Teo [2009] J. M. Karp and G. S. L. Teo, Cell stem cell 4, 206 (2009).
  * Pfeiffer and Höök [2004] I. Pfeiffer and F. Höök, Journal of the American Chemical Society 126, 10224 (2004).
  * Zhang and Winfree [2009] D. Y. Zhang and E. Winfree, Journal of the American Chemical Society 131, 17303 (2009).
  * Markham and Zuker [2005] N. R. Markham and M. Zuker, Nucleic acids research 33, W577 (2005).
  * Markham and Zuker [2008] N. R. Markham and M. Zuker, Bioinformatics: structure, function and applications , 3 (2008).
  * Dimitrov and Zuker [2004] R. A. Dimitrov and M. Zuker, Biophysical Journal 87, 215 (2004).
  * SantaLucia Jr [1998] J. SantaLucia Jr, Proceedings of the National Academy of Sciences 95, 1460 (1998).
  * Jana and Mognetti [2019] P. K. Jana and B. M. Mognetti, Physical Review E 100, 060601 (2019).
  * Gillespie [1977] D. T. Gillespie, The journal of physical chemistry 81, 2340 (1977).
  * Mognetti [2023] B. Mognetti, Model diffusion, https://github.com/bmognetti/Gradient-of-Ligands-/blob/main/ModelDiffusion_LE.ipynb (2023).
  * Martinez-Veracoechea  _et al._ [2014] F. J. Martinez-Veracoechea, B. M. Mognetti, S. Angioletti-Uberti, P. Varilly, D. Frenkel, and J. Dobnikar, Soft matter 10, 3463 (2014).
  * Widder  _et al._ [2022] C. Widder, F. Koch, and T. Schilling, The Journal of Chemical Physics 157 (2022).
  * Netz [2023] R. R. Netz, Derivation of the non-equilibrium generalized langevin equation from a generic time-dependent hamiltonian (2023), arXiv:2310.00748 [cond-mat.stat-mech] .
  * Sforzi  _et al._ [2020] J. Sforzi, L. Palagi, and S. Aime, Biology 9, 202 (2020).
  * Liu and Boyd [2013] Q. Liu and B. J. Boyd, Analyst 138, 391 (2013).
  * Mazur  _et al._ [2017] F. Mazur, M. Bally, B. Städler, and R. Chandrawati, Advances in colloid and interface science 249, 88 (2017).
  * Monteiro  _et al._ [2014] N. Monteiro, A. Martins, R. L. Reis, and N. M. Neves, Journal of the Royal Society Interface 11, 20140459 (2014).
  * Shafiei  _et al._ [2021] M. Shafiei, M. N. M. Ansari, S. I. A. Razak, and M. U. A. Khan, Polymers 13, 2529 (2021).
  * Sarkar and Bose [2019] N. Sarkar and S. Bose, ACS applied materials & interfaces 11, 17184 (2019).
  * Zadeh  _et al._ [2011] J. N. Zadeh, C. D. Steenberg, J. S. Bois, B. R. Wolfe, M. B. Pierce, A. R. Khan, R. M. Dirks, and N. A. Pierce, Journal of computational chemistry 32, 170 (2011).
  * Angelova and Dimitrov [1986] M. I. Angelova and D. S. Dimitrov, Faraday discussions of the Chemical Society 81, 303 (1986).
  * Li [2008] K. Li, The image stabilizer plugin for ImageJ, http://www.cs.cmu.edu/~kangli/code/Image_Stabilizer.html (2008).
  * Fukai and Kawaguchi [2022] Y. T. Fukai and K. Kawaguchi, Bioinformatics , btac799 (2022).
  * Fukai [2021] Y. T. Fukai, laptrack, https://doi.org/10.5281/zenodo.5519537 (2021).