[1]Jennifer L. Ross

[1]Physics Department, Syracuse University, Crouse Drive, Syracuse, 13104, NY,
USA

2]Department of Physics, University of California, Merced, 5200 North Lake Rd,
Merced, 95343, CA, USA

# Effects of Cytoskeletal Network Mesh Size on Cargo Transport

Nimisha Krishnan<EMAIL_ADDRESS>Niranjan Sarpangala
<EMAIL_ADDRESS>Maria Gamez<EMAIL_ADDRESS>Ajay
Gopinathan<EMAIL_ADDRESS><EMAIL_ADDRESS>*  [

###### Abstract

Intracellular transport of cargoes in the cell is essential for the
organization and functioning cells, especially those that are large and
elongated. The cytoskeletal networks inside large cells can be highly complex,
and this cytoskeletal organization can have impacts on the distance and
trajectories of travel. Here, we experimentally created microtubule networks
with varying mesh sizes and examined the ability of kinesin-driven quantum dot
cargoes to traverse the network. Using the experimental data, we deduced
parameters for cargo detachment at intersections and away from intersections,
allowing us to create an analytical theory for the run length as a function of
mesh size. We also used these parameters to perform simulations of cargoes
along paths extracted from the experimental networks. We find excellent
agreement between the trends in run length, displacement, and trajectory
persistence length comparing the experimental and simulated trajectories.

###### keywords:

intracellular transport, microtubule, cytoskeleton, kinesin, cargo transport

##  1 Introduction

The movement and positioning of large objects inside cells requires energy-
using active transport by motor proteins traversing along cytoskeletal
filaments [1]. This process of intracellular transport is responsible for the
organization and reorganization that cells need to survive. Intracellular
transport is especially important in cells that are long and extended, such as
cilia and axons, or particularly crowded and viscous. In mammalian cells,
which are differentiated into a myriad of cell types, diffusion of large
cellular components is impeded by the complex viscoelastic nature of the cell
interior, so active intracellular transport is required.

Cytoskeletal filaments, microtubules and actin, serve as the tracks for
intracellular transport. Microtubules are particularly used for long-distance
transport [1, 2]. Prior works have shown that the arrangement of the
cytoskeletal filaments can affect the transport properties of single motors
and teams of motors attached to cargoes [3, 4, 5, 6, 7, 8, 9, 10, 11]. For
long-distance transport, the microtubules are arranged in logical parallel
bundles to take advantage of kinesin motors that move distally toward the
microtubule plus ends and cytoplasmic dynein motors that move inward to the
microtubule minus ends [12, 13]. In other cell types or locations, the
cytoskeletal networks are more complicated. For instance, in muscle cells,
microtubules create a cross-hatched network creating intersections for
organelles and plasma membrane to anchor during large scale extensions and
contractions [14, 15]. Prior experimental cellular work has demonstrated that
the organization of the cytoskeleton can control the association,
dissociation, and trajectory of vesicles that can dynamically change in time
and space [11, 16].

There are still open questions about how dense, complex, and crowded
conditions can regulate, control, and inhibit intracellular transport. In
order to probe the parameters of control, we created microtubule networks of
varying densities, characterized by the mesh size, the distance between
intersections of microtubules. Using these networks we experimentally probed
the trajectories of kinesin-laden quantum dots as they traverse the network.
The same networks were used as the basis for simulated trajectories for
cargoes where the rates of dissociation at intersections and along the
filaments were determined from experiments, making for closer comparisons. We
also deduced an analytical function for the run length dependence on mesh size
to compare to both experiments and simulations. By comparing the experimental,
simulated, and theoretical results, we determined the effects of the network
mesh size on the motion of kinesin-driven cargoes.

##  2 Methods

###  2.1 Materials and Reagents

Unless otherwise stated, all reagents were purchased from ThermoFisher.

####  2.1.1 Microtubule preparation

Lyophilized 488-tubulin and unlabeled-tubulin were purchased from
Cytoskeleton. Tubulin was resuspended in PEM-80 (80 mM PIPES pH 6.9, 1 mM
MgCl2 and 1 mM EGTA) to a final concentration of 5 mg/ml. Labeled tubulin was
added to unlabeled tubulin at a 1:10 ratio. To polymerize the tubulin, we
added GTP to a final concentration of 1 mM and incubated for 20 min at 34oC.
Finally, we added 20 $\mu$M Taxol to stabilize the polymerized microtubules
and incubated at 34oC for 20 min to equilibrate the Taxol. Microtubules were
kept on the bench and further dilutions of microtubules required 20 $\mu$M
Taxol to keep the filaments stabilized.

####  2.1.2 Kinesin preparation

Kinesin motors were expressed and purified from the pWC2 plasmid available at
AddGene to create a protein with a kinesin-1 motor truncated at amino acid
401, a BCCP tag to allow biotinylation during expression in bacteria, and a
6x-his tag to purify using a nickel affinity column. Kinesin was purified
using standard protocols previously described [17, 18, 19]. Briefly, the
plasmid containing the kinesin construct was transfected into BL21 cells (New
England Biolabs) and bacteria were selected using ampicillin in the media.
Overnight cultures that included biotin in the media were pelleted and the
bacteria was lysed using sonication and chemical lysis. The supernatent was
separated from the bacteria debris using centrifugation and then incubated
with nickel beads to bind the 6x-his tagged protein. Kinesin was eluted using
imidizole and fractions with kinesin were desalted to remove excess imidizole.
Kinesin was aliquoted and snap frozen using liquid nitrogen and stored at
-80oC.

####  2.1.3 Microtubule network preparation

Microtubule networks of varying filament mesh density were made by flow-
aligning microtubules in a crossed-path flow chamber as previously described
[5]. We made the chamber by adhering four square pieces of double-sided tape
on a glass slide such that it made a crossed flow path (Fig. 1A). The slide
was bound to a silanized cover glass treated with hydrophobic silane, PlusOne
Repel silane (Cytiva) as previously described [20].

To create the sample, the following reagents were flowed into the chamber.
First, we flowed 15 $\mu$l of 10$\%$ $\alpha$-tubulin antibody (YL1/2) into
the chamber and incubated for 5 minutes. This surface layer provided a
specific interaction to the microtubules and helped to elevate them above the
polymer surface coating, which was added next. The polymer surface was made by
adding 10 $\mu$l of 5$\%$ Pluronic F-127 block copolymer from both directions
of the flow chamber and incubating for 5 min. The pluronic blocks the surface
from other proteins non-specifically binding. Next we washed the chamber with
wash buffer (90 $\mu$l PEM-80, 10 $\mu$l of 0.5$\%$ Pluronic F-127). Now that
our surface was well coated and blocked, we flowed 10 $\mu$l of polymerized
microtubules diluted to 0.5 mg/ml tubulin concentration. We flowed from the
x-direction, incubated for 2 minutes, washed with wash buffer, and incubated
for another 3 minutes. We repeated the same process in the y-direction. The
chamber was imaged to ensure microtubule networks were bound to the surface
and at the densities needed.

Figure 1:  Microtubule network creation and analysis. (A) Cartoon schematic of
crossed flow path sample chamber and microtubule network. (B) Example image of
microtubule channel for a network created in a crossed channel chamber. (C)
Network image from panel (B) binarized to make a black and white image to be
used for mesh analysis. (D) Skeletonization of network used to automatically
detect intersections and branch lengths. (E) Extracted network used to perform
simulations of motors on networks with the same organization as experiments.
(F) Comparison between the mesh size measured from ImageJ and extracted using
MatLab. Not all networks used in experiments were extracted and used for
simulations. For all images, the scale bar is 5 $\mu$m.

####  2.1.4 Quantum-dot cargo preparation

Quantum dot cargoes were made by mixing streptavidin-labeled quantum dots
(ThermoFisher) with biotinylated kinesin at a ratio of 1:2 and incubated for
one hour on ice (Fig. 2). These cargoes were then diluted by 30 times in
PEM-80 to be used in microtubule networks of varying densities. The final step
of sample chamber preparation was to add the kinesin cargo sample to the
chamber that has been examined on the microscope. The final flow through
contains quantum dots diluted to 1:30 in PEM-80 with 2 mM ATP, 66 mM DTT, and
an oxygen scavenging system, which was 0.66 mg/ml glucose oxidase, 1.5$\%$
final dilution of aqueous catalyse (Sigma catalogue number C30), and 20 mg/ml
glucose in PEM-80. The microtubule network in the crossed part of the flow
chamber had varying densities allowing us to take data in several locations
within the same chamber (Fig. 1). If needed, kinesin cargo sample was
replenished into the same sample chamber to replace the ATP and oxygen
scavenging species that degraded during the assay and allow longer imaging.

Figure 2:  Quantum dot cargo methods and trajectory analysis. (A) Cartoon
schematic of quantum dot cargo attached to a kinesin motor that can walk along
microtubules. (B) Example cargo trajectories (magenta) are displayed on a
dense microtubule network (white) with a small mesh size. Only a subset of
total trajectories are shown. (C) Example cargo trajectories (magenta) along a
sparse microtubule network (white) with a large mesh size. (D) Example
simulated trajectories (magenta) along dense extracted network (white) with a
small mesh size. (E) Example simulated trajectories (magenta) along sparse
extracted network (white) with a large mesh size. For all images, the scale
bar is 5 $\mu$m.

####  2.1.5 Microscopy imaging

Image data was captured with a Nikon Ti-E microscope using epi-fluorescence
and total internal reflection fluorescence (TIRF) microscopy as previously
described [19, 4, 5]. Microtubules were imaged in epi-fluorescence in the
green fluorescence channel using a Hg-Xe illumination source with 480 $\pm$ 25
nm excitation filter, a 500 nm long pass filter, and a 525 $\pm$ 55 nm
emission filter (Chroma). The illumination for the TIRF system was a custom-
built laser system using a 647 nm solid state laser brought into the back of
the 60x, 1.49 NA objective, as previously described [4]. The filter set had no
excitation filter, a 640 nm long pass for the dichroic, and a 680 $\pm$ 50 nm
emission filter (Chroma). All images were made using an IXON electron-
multiplier CCD camera (Andor) with a pixel size of 160 nm. The laser and
camera systems were controlled through Nikon Elements software and images were
recorded to RAM and saved a .nd2 files as uncompressed tif stacks and
metadata. The time series data sets of quantum dot cargoes were taken for 2
mins with 1 s in between frames with an exposure time of 100 ms.

###  2.2 Quantitative image analysis

####  2.2.1 Network mesh size characterization

The control parameter for these studies was the network mesh size, which was
denoted as the distance between neighboring intersections of the microtubule
network. We noticed that the mesh size could change significantly over the
imaging region of our camera, which was 82 $\mu$m on a side. In order to have
the entire region have a similar mesh size, we divided each image into
quarters for analysis of both the mesh size and the trajectories (see below).
This gave more consistent network mesh sizes over an area of 41x41 $\mu
m^{2}$. The same networks were extracted and used as the basis for simulated
trajectories (see below).

We quantified the distance between intersections using the FIJI/ImageJ
AnalyzeSkeleton (2D/3D) plugin [21]. First, we smoothed the images to remove
fluctuations in the background caused by shot noise. Next, we performed
background subtraction on the images to remove global intensity variations due
to the imaging. We converted the image into a binary image by using the auto
threshold function to make the microtubules white on a black background (Fig.
1C). This helps in distinguishing signal (microtubules) from background. We
then skeletonized the image using the binary/skeleton command in FIJI/ImageJ
(Fig. 1D). Finally, we applied the AnalyzeSkeleton (2D/3D) plugin with prune
ends and prune cycle shortest branch enabled. The data reported were the
largest shortest path, detailed info, and labeled skeletons. The resulting
data displayed the branch information, which was saved as a text file. We used
the branch length given in microns as the data for the mesh size. The
statistics of the data were calculated and the distribution was normal, so the
mean and the median were the same. The error bars reported are the standard
error of the mean. The number of branches analyzed are given in Appendix table
1.

####  2.2.2 Tracking and transport analysis

We used Fiji/ImageJ tracking plugin Trackmate [22] to track quantum dot
trajectories. Within Trackmate, we used a setting of 6 pixels for the diameter
of the objects and we allowed a gap in time of 2 frames, so that more than 2
frames without detecting the object nearby resulted in terminating the
measurement. We also used a minimum cut off run length of 3 pixels (160
nm/pixel) and 3 frames (1 s/frame). The localization was allowed to be sub-
pixel. Each of the tracks were manually checked against the movie to ensure
the trajectory tracked was reasonable. The x,y position data over time were
used as the trajectories for further analysis. Example trajectories on
different networks are show in figure 2. The number of tracks analyzed for
each network is given in Appendix table 1.

The run length of a trajectory was determined as the contour length of the
trajectory where the absolute value of all the displacements were summed. The
displacement of a trajectory was defined as the end-to-end distance of the
trajectory. The instantaneous speed was calculated as the positive
displacement between two frames divided by the time between frames. The
average speed was determined as the run length divided by the total time the
quantum dot was associated to the network. For all data types, the data was
averaged and the standard deviation or standard error of the mean was used as
error bars. The number of tracked trajectories used for each network is given
in Appendix table 1.

To quantify the characteristic persistence length of the trajectory, we
calculated the mean squared displacement (MSD) and plotted it versus the
contour length of the trajectory. The MSD was calculated by measuring the
displacement for all points along the trajectory at a specific lag time (time
between frames). For each lag time, the displacement values were squared and
averaged. This is performed for all the lag times that have more than 5 data
points and plotted for only the first 80 points. Once we calculate the MSD as
a function of lag time, we also determined the run length (total contour
length from initial time) of the trajectory for all time. The run length from
zero time is plotted as the x-coordinate and the MSD as a function of lag time
is plotted on the y-coordinate. The MSD was fit to a worm-like chain model
with this equation:

$MSD(L_{c})=2L_{p}L_{c}\left[1-\frac{L_{p}}{L_{c}}\left(1-\exp(-L_{c}/L_{p}\right)\right]$
(1)

where $MSD(L_{c})$ is the mean squared displacement as a function of the
contour length, $L_{c}$, and $L_{p}$ is the persistence length, which is a fit
parameter.

###  2.3 Network extraction and characterization

We used the Matlab tool, FIRE [23] on skeletonized images (Fig. 1D) to extract
corresponding networks (numerical matrices of filaments and vertices) for
cargo transport simulations. We filtered out any filaments that were too short
to be a real filament (typical cut-off length is 0.48 $\mu$ m). We then
characterized the network by quantifying the mesh size, as the mean of the
distances between filament intersections and the persistence length of
filaments. To obtain the persistence length, we restructured the network data
so as to represent filaments as trajectories on a 2D cartesian plane. Then we
measured the MSD of these paths as a function of path length, and fit the
worm-like chain model (Eq. 1) as described above.

###  2.4 Simulations on extracted network

The computational model for cargo transport on the network was similar to
previous works [24, 25]. A cargo of radius $r$ was initialized at a random
point in the 2-dimensional box containing the extracted network. It was
allowed to diffusively search and bind to a filament. The simulation time
starts from the time when the cargo bound to the filament. After binding, the
cargo ballistically moved toward one of the filament ends, i.e. in a time step
$\Delta t$, the cargo position moved by $v\Delta t$ towards the next vertex,
where $v$ is the velocity of cargo. The polarity of the filament was chosen
randomly when the cargo bound and was fixed throughout the given cargo run. As
a cargo walked ballistically on the filament, it could stochastically detach
anywhere along the path with a rate, $k_{off}=v\tilde{k}_{d}$, where
$\tilde{k}_{d}$ was the detachment rate per unit length along a filament. In
addition to stochastic detachments a the cargo walks along a filament, it also
could detach with a fixed probability, $P_{d}$ at filament intersections. A
cargo run stopped when it detached from the filament. Cargo was assumed to
interact with the intersection when it was closer than one cargo radius ($r$)
to one of the filament intersections. The value of the off rate,
$\tilde{k}_{d}$ and detachment probability at intersections, $P_{d}$, were
determined from the analysis of run lengths from manual tracking of
experimental videos (see Sec 3.2). Finally, the cargo was assumed to detach
from filaments when it reached one of the filament ends. In our model, we
simulated only cargo transport along the filament and considered only
ballistic motion for analysis to compare with experimental tracks. This was
different from previous models where diffusive transport of cargo in the
cytoplasm was considered [24, 25] in addition to the ballistic phase on
filaments. Since our goal was to analyze the impact of network features on
transport, which is independent of time, we decided to work with quantities
that do not depend on time. Instead, we examined run length, displacement, and
MSD as a function of path length to compare with experiments. Thus $v$ in our
model was a free parameter that was chosen depending on computational
convenience. We performed N=1000 cargo runs for each network.

##  3 Results and Discussion

###  3.1 Network mesh size

We created artificial cargoes from quantum dots decorated with kinesin-1
motors added to microtubule networks of various mesh sizes. The mesh size of
the networks and a variety of transport parameters were quantified to
determine how the microtubule network density affected the mobility of kinesin
cargoes. The average mesh sizes measured from experimental images using ImageJ
ranged from 1.5 $\mu$m to 5 $\mu$m between intersections (Fig. 1F). The
experimental networks were extracted to be used for simulations. When the mesh
sizes were characterized from simulations using MatLab, they were 37$\%$
higher than the experimental characterization on 97$\%$ of the networks (Fig.
1F). It is possible that in the extraction, we lost small intersections that
were actually there in experiments, resulting in a larger mesh size in
simulations. This could result in altered quantitative results from
simulations compared to experiments. We expected that the trends should be the
same, if the correct underlying mechanisms were being simulated, which we test
below. For plotting of experimental and simulated data, we chose to use the
characteristic mesh sizes determined from experimental images using ImageJ in
all figures.

###  3.2 Run length

The run length of a quantum dot cargo is the total distance, or contour
length, traveled along the microtubules of the network before detaching.
Experimental data was collected and long trajectories were manually tracked,
regardless of network mesh size. The run length, $s$, was determined from the
manually tracked trajectories and the data was binned with either 8 $\mu$m or
1 $\mu$m bins (Fig. 3A, blue). We deduced the two parameters that controlled
run length via dissociation, namely the off rate for cargoes between
intersections, $k_{off}$, and the probability of detaching at the
intersections, $P_{d}$, by fitting the run length histograms to exponential
decay functions of the form: $y=ae^{-(\tilde{k_{d}}+P_{d}/\lambda)s}$, where
$\lambda$ is the mean mesh size for all extracted networks, $a$ is an
arbitrary normalization parameter, and $\tilde{k_{d}}$ is the detachment rate
per unit length. The parameter $\tilde{k_{d}}$ is equivalent to the inverse of
the “natural” run length for the cargoes in the absence of intersections and
is equal to the off rate between intersections divided by the cargo velocity:
$k_{off}/v$ (Fig. 3A, blue).

Figure 3:  Determination of simulation parameters and analytical theory. (A)
Distribution of run lengths of all tracks (blue circles, N=246) fit to
exponential decay (blue line). Distribution of run lengths of tracks that did
not visit an intersection (red circles, N=20) fit to exponential decay (red
line). Short run lengths were excluded from the fit because manual tracking
had a systematic bias against short run lengths. (B) Analytical theory for the
average run length as a function of mesh size (Eq. 2) for various values of
$P_{d}$. The legend indicates the value of $P_{d}$ plotted with the two values
from the histograms, 0.62 and 0.75 in gray.

Because we have two unknown parameters, we also need to use a second set of
data to deduce the off rate between intersections. Using only the subset of
tracks that never visited an intersection during their trajectories, the run
lengths were again binned and fit to $y=ae^{-(\tilde{k_{d}}+1/\lambda)s}$
(Fig. 3A, red). The two histograms had fit parameters of
$A=\tilde{k_{d}}+P_{d}/\lambda$ and $B=\tilde{k_{d}}+1/\lambda$, which were
used to deduce $\tilde{k_{d}}$ and $P_{d}$ (Table 4).

The histogram bin size had an effect on the fitting parameters and hence the
values of $k_{off}$ and $P_{d}$ used for simulations. Histograms with 8 $\mu$m
bins resulted in a $P_{d}$ value of 0.75 and histograms with 1 $\mu$m bins
resulted in a $P_{d}$ value of 0.62. Assuming the distribution of run lengths
has this form: $ae^{-(\tilde{k_{d}}+P_{d}/\lambda)s}$, then the average run
length, $\langle s\rangle$, for a network with a mesh size, $\lambda$, should
be given by finding the average of this expression given by:

$\langle s\rangle=1/(\tilde{k_{d}}+P_{d}/\lambda).$ (2)

Thus the dependence of the average run length ($\langle s\rangle$) with mesh
size ($\lambda$) is non-linear and saturates to $1/\tilde{k_{d}}$ for large
mesh sizes (Fig. 3B). Given the uncertainty of the estimation of $P_{d}$, we
can examine the sensitivity of this analytical expression to the value of
$P_{d}$ (Fig. 3B). For small $P_{d}$, the run length saturates to the natural
run length for cargoes on a single microtubule. When $P_{d}$ approaches 1, the
run lengths are depressed and only reach the natural run length at higher mesh
sizes (Fig. 3B). This theoretical equation is general in that it allows us to
compute run lengths at any values of motor off rate, mesh size, and detachment
probability at intersections. This equation will also allow future experiments
to compute the detachment probability at intersections without having to
manually filter trajectories that encounter intersections. For the
simulations, we used a value of $P_{d}$ equal to 0.75 and a $\tilde{k_{d}}$
equal to 0.1.

We compared our analytical results to both our experimental quantification of
run length and simulations of the run length on our microtubule networks. For
the experimental data, we observed that as the mesh size increased, the
average run length also increased (Fig. 4Ai,ii). We assigned an artificial
cut-off between low mesh size (less than 2 $\mu m$) and high mesh size (larger
than 2 $\mu m$). Using this cut-off, we found that the median run length for
low mesh size was 3.0 $\pm$ 0.1 $\mu m/s$ and the high mesh size was 3.9 $\pm$
0.1 $\mu m/s$ (Fig. 4Aii).

Figure 4:  Total run lengths for experimental and simulated trajectories. (A)
Experimental results. (i) Plot showing the average run length ($\mu$m) for
each network of a given average mesh size ($\mu$m). Best fit slope is 0.5
$\pm$ 0.1. (ii) Comparison of the distribution of the average run lengths for
networks with low mesh size less and high mesh size. (B) Simulation results.
(i) Average run lengths of simulated trajectories on the same networks. Each
data point represents average over N=1000 cargo runs, error bars represent
standard error of the mean. Best fit slope is 0.4 $\pm$ 0.1. (ii) Comparison
of the distribution of the average run lengths for low mesh size and high mesh
size. All fit parameters for data are given in Appendix table 2. (C)
Comparison of experimental run lengths (magenta circles) and simulated run
lengths (blue circles) as a function of network mesh size with the analytical
theory with various values of $P_{d}$, given in the legend.

There was a distinct difference in the statistics for networks with mesh sizes
above and below this cut-off (Fig. 4Aii). Indeed, performing the Kolmogorov-
Smirnoff statistical test (KS Test), we found the probability was p = 0.0003
or 0.03% that the small and large mesh size results were the same (Fig. 4Aii).
Thus, we concluded that the threshold at a mesh size of 2 $\mu m$ is a
reasonable cut-off between low and high mesh sizes for further comparisons.
The experimental run length results implied that quantum dots cover longer
distances when the microtubule tracks were more open, with fewer
intersections. This result makes sense, since smaller mesh sizes should have
more intersections, and kinesin has been shown to have a high probability of
terminating a run when contacting a microtubule intersection [5].

Using the deduced probabilities for dissociation between or at intersections
(Fig. 3A), we were able to simulate trajectories on different networks and
quantify the run lengths as a function of mesh size (Fig. 4B). The simulated
run lengths increased with increasing mesh size just like the experimental
data with a similar slope (Fig. 4Bi). Using the same mesh size cut-off, we
compared the simulated run lengths to find the low mesh size median run length
was 4.2 $\pm$ 0.8 and the high mesh size median run length was 4.8 $\pm$ 0.1.
These results were significantly different with a probability of 0.0015% that
they are the same distribution using the KS Test.

We can compare the experimental and simulated trajectory run lengths to the
analytical expression for average run length as a function of mesh size (Eq.
2, Fig. 4C). Plotting all together, it is clear that the simulations have
systematically longer run lengths than the experimental results. Given the
sensitivity to the value of $P_{d}$, it is possible that adjusting this
parameter or the natural run length, $1/\tilde{k_{d}}$, could cause the
difference. Examining the data with multiple values of $P_{d}$ plotted, we
estimate that the probability of detaching at an intersection is between 0.3
and 0.75 for the experimental data and between 0.2 and 0.5 for the simulation
data (Fig. 4C). This is a not surprising considering that the extracted
networks have larger mesh sizes compared to experiments (Fig. 1). Thus, for
the same values of $P_{d}$ and $1/\tilde{k_{d}}$, the number of intersections
encountered is smaller, and the run lengths will thus be larger.

Another difference between simulations and experiment could results from track
switching at intersections. In our simulations, we did not include track
switching as an option at intersections (Table 4) because switching has been
shown to be infrequent for kinesin cargoes with one or two motors [5].
Further, in our manual tracking of long trajectories, we only observed
switching at intersections with about 5% probability, which matches prior
reports for single GFP-kinesin [5]. In order to check if 5% switching
probability at intersections could alter the results, we simulated
trajectories in the networks including this probability. We found that this
small switching probability has no effect on the run lengths we observe in
simulations (Appendix Fig. 81), justifying our choice to not include it in the
simulations. Further, this also implies that there is likely on average one
active motor per quantum dot, although there is a small probability that some
quantum dots have two motors.

###  3.3 Instantaneous and average speeds

We used metrics that were independent of time to allow us to compare between
the experiment and simulation data. In our model, the speed $v$ was a free
parameter that was chosen depending on computational convenience. Given these
assumptions, we would expect that the average and instantaneous velocities
should be independent of mesh size. As a rationale check, we calculated the
instantaneous speed of the quantum dot cargoes between two time points. We
found that there was a shallow trend in the median of the instantaneous speed
as a function of average mesh size (Fig. 5A). Using the same threshold at a
mesh size of 2 $\mu m$, we found that the two distributions have different
medians of 0.052 $\pm$ 0.004 $\mu m/s$ for the low mesh size and 0.060 $\pm$
0.002 $\mu m/s$ for the high mesh size (Fig. 5B). The distribution of
instantaneous speeds was wider with a standard deviation of 0.018 $\pm$ 0.004
compared to the standard deviation of the high mesh size data, which had a
standard deviation 0.010 $\pm$ 0.002 (Fig. 5B). This was likely because, at a
smaller mesh size, there was a higher probability of a cargo encountering an
intersection between two frames and affecting the instantaneous speed.
Comparing the data with the KS Test, we found that the distribution in
instantaneous velocities for the low and high mesh sizes were not
statistically different, with a 5.6% probability that they are the same.

We also quantified the average velocity of the cargo trajectories, given by
the total run length (Fig. 4) over the total association time. Since the
network intersections reduced both the association time and the run length by
the same mechanism, specifically causing the dissociation of the cargoes, we
expected both parameters to decrease similarly. Indeed, we found that the
average speed of the cargoes was unaffected by mesh size, as expected (Fig.
5C). The average speed distributions had medians of 0.82 $\pm$ 0.06 for low
mesh and 0.80 $\pm$ 0.03 for high mesh size (Fig. 5D). The probability that
these two distributions are the same is 39% using the KS test (Fig. 5D). These
checks indicated that our model assumptions about velocity were reasonable.

Figure 5:  Cargo speeds were independent of mesh size. (A) Plot showing the
median instantaneous speed ($\mu$m/s) for each network of a given average mesh
size ($\mu$m). The best fit slope was 0.005 $\pm$ 0.003. (B) Comparison of the
distribution of the median instantaneous speed for low and high mesh sizes,
which were statistically the same. (C) Plot showing the median average speed
($\mu$m/s) for each network of a given average mesh size ($\mu$m) with the
best fit slope of -0.03 $\pm$ 0.03. (D) Comparison of the distribution of the
median average speed for low and high mesh sizes, which were statistically the
same. All fit parameters for data are given in Appendix table 2.

###  3.4 Displacement and tortuosity

The displacement is the end-to-end length of the cargo’s trajectory. We
quantified the displacement for experimentally measured trajectories and found
that it increased with the mesh size (Fig. 6Ai). The median displacement for
low mesh size was 1.77 $\pm$ 0.07 and for high mesh size was 2.26 $\pm$ 0.08
(Fig. 6Aii). The probability that the distributions in displacement were the
same was 0.01% using the KS Test.

The simulated trajectories also showed the same trend in displacement,
increasing with mesh size (Fig. 6Bi). The median displacement for low mesh
size was 3.12 $\pm$ 0.06 and for high mesh size was 3.80 $\pm$ 0.01 (Fig.
6Bii). The probability that the distributions in displacement were the same
was 0.0015% for high and low mesh sizes using the KS Test. Although the trends
are similar, the absolute numbers for the simulated displacements were higher
than the measured displacements for the same networks. Like the run length,
the displacements of simulated trajectories were not affected by a 5%
probability to switch microtubules at intersections (Appendix Fig. 81).

Figure 6:  Total displacement for experimental and simulated trajectories. (A)
Displacement of experimental trajectories. (i) Plot showing the average
displacement ($\mu$m) against average mesh size ($\mu$m). The best fit slope
is 0.35 $\pm$ 0.06. (ii) Comparison of the distribution of average
displacements for low mesh size and high mesh size. (B) Displacement of
simulated trajectories. (i) Plot showing the average displacement ($\mu$m)
with average mesh size ($\mu$m). The best-fit slope is 0.5 $\pm$ 0.1. (ii)
Comparison of the distribution of average displacements for low mesh size and
high mesh size. All fit parameters for data are given in Appendix table 2.

We can determine the tortuosity of the trajectories by dividing the contour
length by the displacement. This is a parameter used for examining the flow of
material through porous media, and can be used to characterize the mobility.
We find that the average tortuosity for low mesh size networks is 1.7 $\pm$
0.1 and for high mesh size is 1.6 $\pm$ 0.1, which are the same. Thus, the run
length and displacement are being rescaled within the network by the same
process, most likely the presence of the intersections.

####  3.4.1 Mean square displacement

Figure 7:  Mean squared displacement and persistence length of trajectories.
(A) Mean squared displacement for experimental data (i) All experimental data
sets from low mesh size (cyan circles) and high mesh size (green circles)
plotted together. These data were fit to the worm-like chain model (Eq. 1) to
find the persistence length. (ii) Persistence length ($\mu$m) plotted against
mesh size ($\mu$m). The best fit slope is 0.2 $\pm$ 0.1 (iii) Comparison of
median of persistence lengths for low mesh size and high mesh size. (B) Mean
squared displacement for simulation data (i) All data sets from low mesh size
(cyan circles) and high mesh size (green circles) plotted together. These data
were fit to the worm-like chain model (Eq. 1) to find the persistence length.
(ii) Persistence length ($\mu$m) plotted against mesh size ($\mu$m). The best
fit slope is 1.6 $\pm$ 0.2 (iii) Comparison of average of persistence length
for low mesh size and high mesh size simulations. All fit parameters can be
found in Appendix table 2.

The characteristic persistence length of the trajectories of the quantum dot
cargoes was characterized using a mean squared displacement (MSD) as a
function of trajectory contour length. Fitting the data this way removed the
need to consider time dependence, as described above. Each MSD was fit to a
worm-like chain model (Eq. 1) to determine the persistence length, $L_{p}$ of
the trajectory (Fig. 7Ai,Bi). For both experimental data and simulated data,
the MSD was cut off at a countour length of 6 $\mu$m for the fitting in order
to compare them.

For the experimental trajectories, the persistence length depended on the mesh
size linearly (Fig. 7Aii). The median for low mesh size was 0.44 $\pm$ 0.06
and the median for high mesh size was 0.8 $\pm$ 0.2. The standard deviation
for the high mesh size data was large (SD = 0.7 $\pm$ 0.2) compared to the low
mesh size data (SD = 0.25 $\pm$ 0.06), which resulted in a probability of
4.2$\%$ that these distributions are the same using the KS Test (Fig. 7Aiii).

The simulation trajectory data showed the same linear trend with mesh size
(Fig. 7Bii), but the data was less spread out for both small and large mesh
sizes. The median for low mesh size was 2.8 $\pm$ 0.3 and the median for high
mesh size was 4.1 $\pm$ 0.3. The KS test revealed that the probability they
are the same distribution is only 0.3$\%$ (Fig. 7Biii).

##  4 Conclusion

The organization of cytoskeletal filaments is known to affect the motion of
motors and cargoes. Here, we took the approach to characterize the network
using the mesh size and examine the effects of the mesh size on the run
length, displacement, and mean square displacement of the motion through the
network. We found that these parameters are sensitive to mesh size, even over
the small range in mesh size that we are able to realize in experiments, 1
$\mu$m to 5 $\mu$m. Despite the small range of mesh sizes, there were
significant changes in all the trajectory parameters.

Using experimental trajectories, we deduced the parameters for cargo
detachment at intersections and between intersections. We used these
parameters to create an analytical theory for the average run length as a
function of mesh size, which had similar trends as our experimental data and
was sensitive to the probability of detaching at intersections.

Using the exact same networks extracted from the experimental data and the off
rates, we were able to simulate cargo trajectories through the networks. The
simulation data had the same trends and similar quantitative results as the
experiments for run length, displacement, and mean squared displacement for
real cargoes assembled with kinesin motors. We anticipate that future work
with different mesh sizes, filament organizations, and motor types will be
modeled with the same fundamental principles we uncover here. Specifically, we
anticipate that other motors would have different reactions to intersections
and composite motor systems would further increase complexity.

Acknowledgments We would like to acknowledge the support of the members of the
Ross Lab in the Physics Department at Syracuse University. NK was partially
supported on funds from Syracuse University, the National Science Foundation
grant NSF BIO-2134215 to JLR and National Science Foundation grant
DMREF-2118403 to JLR. AG, MG and NS acknowledge support from the National
Science Foundation (NSF-DMS-1616926 to AG) and NSF-CREST: Center for Cellular
and Biomolecular Machines at UC Merced (NSF-HRD-1547848 and NSF-HRD-2112675 to
AG). AG and NS also acknowledge support from the NSF Center for Engineering
Mechanobiology grant CMMI-1548571. NS acknowledges Graduate Dean’s
dissertation fellowship from the University of California, Merced.

## Declarations

  * •

Funding: The research leading to these results received funding from Syracuse
University and the National Science Foundation grants NSF BIO-2134215 and
DMREF-2118403 as well as NSF-DMS-1616926, NSF-HRD-1547848, NSF-HRD-2112675,
NSF-CMMI-1548571 and UC Merced Graduate Division. Employment: JLR and NK are
employed by Syracuse University. NS and AG are employed by University of
California, Merced.

  * •

Data Availability: Microscopy and simulation data generated for this
manuscript are available upon request from the corresponding authors JLR for
experimental data and AG for simulation data.

  * •

Authors’ contributions: NK conceived and design the work, performed data or
the acquisition, analysis, and interpretation of data, drafted and edited the
manuscript, and is accountable for the work. NS performed simulations and
theoretical calculations, analyzed and interpreted data, drafted and edited
the manuscript and is accountable for the work. MG performed simulations and
analyzed and interpreted data. AG performed theoretical calculations,
interpreted data, drafted and edited the manuscript and is accountable for the
work. JLR designed the experimental work, analyzed and interpreted data,
drafted and edited the manuscript and is accountable for the work.

##  Appendix A Experimental Data Appendix

In this appendix, we list the parameters of the data for the trajectories for
each experimental movie.

Table 1:  Summary of mesh size data from experimental networks Network name |  Number of branches |  Average mesh size |  Mesh size standard Error |  Number of tracked trajectories
---|---|---|---|---
HDN001.1 |  566 |  1.675 |  0.025 |  21
HDN001.2 |  620 |  1.73 |  0.029 |  41
HDN001.3 |  423 |  1.839 |  0.039 |  61
HDN001.4 |  601 |  1.676 |  0.024 |  50
HDN002.1 |  581 |  1.678 |  0.026 |  20
HDN002.2 |  536 |  1.784 |  0.035 |  19
HDN002.3 |  418 |  1.882 |  0.039 |  40
HDN002.4 |  469 |  1.765 |  0.033 |  36
HDN005.1 |  1001 |  1.96 |  0.032 |  120
HDN005.2 |  1213 |  1.928 |  0.028 |  115
HDN005.3 |  994 |  1.687 |  0.02 |  185
HDN008.1 |  148 |  2.537 |  0.142 |  24
HDN008.2 |  89 |  3.816 |  0.35 |  40
HDN008.3 |  131 |  3.021 |  0.178 |  56
LDN001.2 |  181 |  2.029 |  0.084 |  34
LDN001.3 |  138 |  2.554 |  0.19 |  32
LDN001.4 |  209 |  1.884 |  0.065 |  41
LDN002.1 |  106 |  2.67 |  0.202 |  41
LDN002.2 |  103 |  2.783 |  0.21 |  53
LDN002.3 |  182 |  2.675 |  0.182 |  44
LDN002.4 |  136 |  2.372 |  0.19 |  49
LDN004.1 |  107 |  3.049 |  0.255 |  69
LDN004.2 |  203 |  2.188 |  0.098 |  71
LDN004.3 |  114 |  3.764 |  0.293 |  117
LDN004.4 |  183 |  2.218 |  0.153 |  108
LDN005.1 |  76 |  3.643 |  0.296 |  36
LDN005.2 |  122 |  3.309 |  0.235 |  41
LDN005.3 |  136 |  3.694 |  0.277 |  56
LDN005.4 |  99 |  3.076 |  0.269 |  42
LDN007.1 |  90 |  3.496 |  0.288 |  48
LDN007.2 |  48 |  4.002 |  0.602 |  28
LDN007.3 |  91 |  3.146 |  0.284 |  52
LDN007.4 |  71 |  4.065 |  0.559 |  63
LDN008.1 |  75 |  2.202 |  0.177 |  12
LDN008.3 |  44 |  4.517 |  0.71 |  37
VHD001.1 |  924 |  1.663 |  0.019 |  82
VHD001.2 |  1003 |  1.687 |  0.019 |  73
VHD001.3 |  858 |  1.696 |  0.022 |  177
VHD002.4 |  1067 |  1.673 |  0.019 |  45
VHD003.3 |  718 |  1.663 |  0.022 |  64
VHD006.3 |  700 |  1.602 |  0.019 |  37
Table 2:  The fit parameters for the linear fits to the data as a function of mesh size for all figures, as denoted.  Y-axis parameter |  Intercept |  Slope |  R-squared |  Chi-squared |  Reference to Figure
---|---|---|---|---|---
Run length - Experiments |  2.19$\pm$0.27 |  0.52$\pm$0.1 |  0.4 |  11.81 |  Fig. 3A(i)
Run length - Simulations |  3.42$\pm$0.26 |  0.42$\pm$0.1 |  0.38 |  6.31 |  Fig. 3B(i)
Instantaneous speed |  0.05$\pm$0.01 |  0$\pm$0 |  0.08 |  0.01 |  Fig. 5A
Average speed |  0.88$\pm$0.09 |  -0.03$\pm$0.03 |  0.02 |  1.33 |  Fig. 5C
Average displacement - Experiments |  1.21$\pm$0.17 |  0.35$\pm$0.06 |  0.44 |  4.64 |  Fig. 6A(i)
Average displacement - Simulations |  2.17$\pm$0.21 |  0.54$\pm$0.08 |  0.60 |  4.1 |  Fig. 6B(i)
Persistence length - Experiments |  0.23$\pm$0.3 |  0.22$\pm$0.12 |  0.09 |  13.12 |  Fig. 6A(i)
Persistence length - Simulations |  -2.69$\pm$2.51 |  1.56 $\pm$ 0.19 |  0.22 |  58.9 |  Fig. 6B(i)

##  Appendix B Simulation Appendix

Table 3:  Summary of each mesh size data from simulation results Network name |  Number of branches |  Average mesh size |  Mesh size standard Error |  Number of tracked trajectories
---|---|---|---|---
H1.1 |  152 |  2.33 |  0.15 |  1000
H1.2 |  154 |  2.40 |  0.14 |  1000
H1.3 |  199 |  2.09 |  0.11 |  1000
H1.4 |  149 |  2.36 |  0.15 |  1000
H2.1 |  235 |  2.20 |  0.09 |  1000
H2.2 |  178 |  2.39 |  0.14 |  1000
H2.3 |  228 |  2.36 |  0.10 |  1000
H5.1 |  230 |  2.34 |  0.13 |  1000
H5.3 |  292 |  2.23 |  0.09 |  1000
H8.1 |  65 |  3.32 |  0.36 |  1000
H8.2 |  60 |  3.33 |  0.38 |  1000
H8.3 |  74 |  3.31 |  0.33 |  1000
H8.4 |  45 |  4.38 |  0.71 |  1000
L1.2 |  27 |  5.58 |  1.10 |  1000
L1.3 |  18 |  4.00 |  1.26 |  1000
L1.4 |  16 |  5.14 |  1.24 |  1000
L2.1 |  43 |  3.21 |  0.28 |  1000
L2.2 |  38 |  3.16 |  0.52 |  1000
L2.3 |  43 |  3.89 |  0.42 |  1000
L2.4 |  13 |  5.20 |  1.32 |  1000
L4.1 |  34 |  4.43 |  0.75 |  1000
L4.2 |  41 |  2.99 |  0.40 |  1000
L4.3 |  54 |  3.84 |  0.51 |  1000
L4.4 |  42 |  4.68 |  1.04 |  1000
L5.1 |  44 |  3.91 |  0.54 |  1000
L5.2 |  25 |  4.57 |  0.76 |  1000
L5.3 |  46 |  4.31 |  0.49 |  1000
L5.4 |  24 |  5.26 |  1.15 |  1000
L7.1 |  34 |  3.61 |  0.60 |  1000
L7.2 |  2 |  7.63 |  6.09 |  1000
L7.3 |  41 |  3.86 |  0.57 |  1000
L7.4 |  27 |  5.15 |  1.27 |  1000
V1.1 |  322 |  2.10 |  0.07 |  1000
V1.2 |  252 |  2.18 |  0.10 |  1000
V1.3 |  266 |  2.09 |  0.08 |  1000
V1.4 |  252 |  2.25 |  0.11 |  1000
Table 4:  Parameters used for running the simulation Description |  Units |  Variable Name |  Value
---|---|---|---
Number of cargoes |  - |  numCargos |  1000
Velocity |  ($\mu m$/s) |  v |  0.1
Cargo radius |  $\mu m$ |  cRad |  0.1
Cargo step size |  ($\mu m$) |  dstep |  0.2
Maximum measuring time |  (s) |  tMax |  1000
Diffusion constant |  ($\mu m$2/s) |  D |  1
Rate of detachment from filament |  (s-1) |  $k_{off}$ |  0.01
Rebinding allowed (1 if allowed, 0 if not allowed) |  - |  reb_all |  0
Probability of detachment at intersections |  - |  $P_{d}$ |  0.754
Inverse of the run length without intersections |  $(\mu m^{-1})$ |  $\tilde{k_{d}}$ |  0.1
Figure 8: Comparison between simulations with different probabilities of
switching at intersections ($P_{s}$). We computed the mean of run lengths and
displacements of simulated cargo runs over all extracted networks. No
statistically significant difference was seen in mean run length or
displacement values between $P_{s}=0$ and $P_{s}=0.05$. Statistical
significance was determined using students t-test, conventions used: n.s. for
$p>0.05$, **** for p$<10^{-4}$.

## References

  *   * Hirokawa [1998] Hirokawa, N.: Kinesin and dynein superfamily proteins and the mechanism of organelle transport. Science 279(5350), 519–526 (1998)
  * Vale [2003] Vale, R.D.: The molecular motor toolbox for intracellular transport. Cell 112(4), 467–480 (2003)
  * Ross et al. [2008] Ross, J.L., Ali, M.Y., Warshaw, D.M.: Cargo transport: molecular motors navigate a complex cytoskeleton. Current opinion in cell biology 20(1), 41–47 (2008)
  * Ross and Dixit [2010] Ross, J.L., Dixit, R.: Multiple color single molecule tirf imaging and tracking of maps and motors. Methods in cell biology 95, 521–542 (2010)
  * Ross et al. [2008] Ross, J.L., Shuman, H., Holzbaur, E.L., Goldman, Y.E.: Kinesin and dynein-dynactin at intersecting microtubules: motor density affects dynein function. Biophysical journal 94(8), 3115–3125 (2008)
  * Ali et al. [2016] Ali, M.Y., Vilfan, A., Trybus, K.M., Warshaw, D.M.: Cargo transport by two coupled myosin va motors on actin filaments and bundles. Biophysical journal 111(10), 2228–2240 (2016)
  * Lombardo et al. [2017] Lombardo, A.T., Nelson, S.R., Ali, M.Y., Kennedy, G.G., Trybus, K.M., Walcott, S., Warshaw, D.M.: Myosin va molecular motors manoeuvre liposome cargo through suspended actin filament intersections in vitro. Nature communications 8(1), 15692 (2017)
  * Lombardo et al. [2019] Lombardo, A.T., Nelson, S.R., Kennedy, G.G., Trybus, K.M., Walcott, S., Warshaw, D.M.: Myosin va transport of liposomes in three-dimensional actin networks is modulated by actin filament density, position, and polarity. Proceedings of the National Academy of Sciences 116(17), 8326–8335 (2019)
  * Walcott and Warshaw [2022] Walcott, S., Warshaw, D.: Modeling myosin va liposome transport through actin filament networks reveals a percolation threshold that modulates transport properties. Molecular Biology of the Cell 33(2), 18 (2022)
  * Nelson et al. [2009] Nelson, S.R., Ali, M.Y., Trybus, K.M., Warshaw, D.M.: Random walk of processive, quantum dot-labeled myosin va molecules within the actin cortex of cos-7 cells. Biophysical journal 97(2), 509–518 (2009)
  * Tabei et al. [2013] Tabei, S.A., Burov, S., Kim, H.Y., Kuznetsov, A., Huynh, T., Jureller, J., Philipson, L.H., Dinner, A.R., Scherer, N.F.: Intracellular transport of insulin granules is a subordinated random walk. Proceedings of the National Academy of Sciences 110(13), 4911–4916 (2013)
  * Verhey et al. [2011] Verhey, K.J., Kaul, N., Soppina, V.: Kinesin assembly and movement in cells. Annual review of biophysics 40, 267–288 (2011)
  * Burute and Kapitein [2019] Burute, M., Kapitein, L.C.: Cellular logistics: unraveling the interplay between microtubule organization and intracellular transport. Annual review of cell and developmental biology 35, 29–54 (2019)
  * Robison et al. [2016] Robison, P., Caporizzo, M.A., Ahmadzadeh, H., Bogush, A.I., Chen, C.Y., Margulies, K.B., Shenoy, V.B., Prosser, B.L.: Detyrosinated microtubules buckle and bear load in contracting cardiomyocytes. Science 352(6284), 0659 (2016)
  * Uchida et al. [2022] Uchida, K., Scarborough, E.A., Prosser, B.L.: Cardiomyocyte microtubules: control of mechanics, transport, and remodeling. Annual Review of Physiology 84, 257–283 (2022)
  * Conway et al. [2014] Conway, L., Gramlich, M.W., Ali Tabei, S., Ross, J.L.: Microtubule orientation and spacing within bundles is critical for long-range kinesin-1 motility. Cytoskeleton 71(11), 595–610 (2014)
  * Gilbert and Johnson [1993] Gilbert, S.P., Johnson, K.A.: Expression, purification, and characterization of the drosophila kinesin motor domain produced in escherichia coli. Biochemistry 32(17), 4677–4684 (1993)
  * Pierce and Vale [1998] Pierce, D.W., Vale, R.D.: [14] assaying processive movement of kinesin by fluorescence microscopy. In: Methods in Enzymology vol. 298, pp. 154–171. Elsevier, ??? (1998)
  * Conway et al. [2012] Conway, L., Wood, D., Tüzel, E., Ross, J.L.: Motor transport of self-assembled cargos in crowded environments. Proceedings of the National Academy of Sciences 109(51), 20814–20819 (2012)
  * Dixit and Ross [2010] Dixit, R., Ross, J.L.: Studying plus-end tracking at single molecule resolution using tirf microscopy. Methods in cell biology 95, 543–554 (2010)
  * Lee et al. [1994] Lee, T.-C., Kashyap, R.L., Chu, C.-N.: Building skeleton models via 3-d medial surface axis thinning algorithms. CVGIP: graphical models and image processing 56(6), 462–478 (1994)
  * Tinevez et al. [2017] Tinevez, J.-Y., Perry, N., Schindelin, J., Hoopes, G.M., Reynolds, G.D., Laplantine, E., Bednarek, S.Y., Shorte, S.L., Eliceiri, K.W.: Trackmate: An open and extensible platform for single-particle tracking. Methods 115, 80–90 (2017)
  * STEIN et al. [2008] STEIN, A.M., VADER, D.A., JAWERTH, L.M., WEITZ, D.A., SANDER, L.M.: An algorithm for extracting the network geometry of three-dimensional collagen gels. Journal of Microscopy 232(3), 463–475 (2008) https://doi.org/10.1111/j.1365-2818.2008.02141.x https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1365-2818.2008.02141.x
  * Ando et al. [2015] Ando, D., Korabel, N., Huang, K.C., Gopinathan, A.: Cytoskeletal network morphology regulates intracellular transport dynamics. Biophysical Journal 109(8), 1574–1582 (2015) https://doi.org/10.1016/j.bpj.2015.08.034
  * Maelfeyt et al. [2019] Maelfeyt, B., Tabei, S.M.A., Gopinathan, A.: Anomalous intracellular transport phases depend on cytoskeletal network features. Phys. Rev. E 99, 062404 (2019) https://doi.org/10.1103/PhysRevE.99.062404


# Spin-orbit coupling induced by ascorbic acid crystals

Florence Grenapin  Department of physics, University of Ottawa, Advanced
Research Complex, 25 Templeton Street, K1N 6N5, Ottawa, ON, Canada  Alessio
D’Errico<EMAIL_ADDRESS>Department of physics, University of Ottawa,
Advanced Research Complex, 25 Templeton Street, K1N 6N5, Ottawa, ON, Canada
Ebrahim Karimi  Department of physics, University of Ottawa, Advanced Research
Complex, 25 Templeton Street, K1N 6N5, Ottawa, ON, Canada  National Research
Council of Canada, 100 Sussex Drive, K1A 0R6, Ottawa, ON, Canada

(June 2022)

###### Abstract

Some anisotropic materials form semicristalline structures, called
spherulites, which observed in a polarisation microscope, exhibit a
characteristic “maltese-cross”-like pattern. While this observation has been
hitherto considered as a tool to characterize these materials, we show that
these patterns are associated with a strong light’s spin-orbit coupling
induced by the spherulite structures. We experimentally demonstrate these
effects using samples of crystallized ascorbic acid and observing the creation
of optical vortices in transmitted laser beams, as well as the formation of
inhomogeneous polarisation patterns. Our findings suggest the use of
spherulites in frequency ranges, e.g. in the THz domain, where polarisation
and spatial shaping of electromagnetic radiation is still a challenging task.

##  I Introduction

Polymer spherulites are spherically symmetric semicristalline structures
typically observed when a molten polymer is slowly cooled down Crist and
Schultz (2016). The slow cooling allows polymer chains to form in ordered
configurations. The crystallization starts around point defects and consists
in lamellae structures that, in the absence of temperature gradients, grow
radially from the defect centre.
The radially directed fibrillar structure can be easily seen at low
magnification. Such fibrils are composed of one or multiple crystals elongated
in the radial direction. As a consequence of the optical anisotropy of the
crystals, spherulites can be considered as one of the simplest examples of
patterned anisotropic media. Such media have been subject of intense research
in the last decade due to their action on optical beams that couples optical
polarisation and inhomogeneous phase transformations. More specifically,
Pancharatnam-Berry Optical Elements (PBOEs) Marrucci  _et al._ (2006a) are
slabs of patterned anisotropic materials whose inhomogeneous orientation of
the optic axis can be used to structure the phase and polarisation of light in
a way that is conditioned by the polarisation of the input beam. The most
known and widely used device of this kind is the $q$-plate, a PBOE with an
optic axis that varies linearly with the azimuthal angle Marrucci  _et al._
(2006b). $Q$-plates have been used as a tool in different areas Rubano  _et
al._ (2019): from quantum information and simulation Nagali  _et al._ (2009);
Cardano  _et al._ (2015) to microscopy and surface structuring Yan  _et al._
(2015); Nivas  _et al._ (2015). Typically $q$-plates are based on liquid
crystals, in which case their retardation can be finely controlled either
thermally Karimi  _et al._ (2009) or electrically Piccirillo  _et al._ (2010).
More recently $q$-plates, as well as other PBOEs, based on subwavelength
structures have been developed: from nanostructured glasses Beresna  _et al._
(2011) to plasmonic metasurfaces Karimi  _et al._ (2014).
Here, we show how slabs of spherulites – specifically how thin layers of
ascorbic acid spherulites – also behave like PBOEs with an azimuthal optic
axis dependence. Indeed, when observed in the context of polarized light
microscopy, i.e. between crossed polarizers, spherulites display “maltese-
cross”-like patterns which strongly resemble $q$-plate patterns with
topological charge $q=1$ (see Fig. 1).
In $q$-plates, the radial optic pattern introduces an azimuthal phase factor
to the light transmitted through it, a factor which is associated with a well-
defined amount of Orbital Angular Momentum (OAM). This process is known as
optical spin to orbit coupling. Depending on the input polarisation state one
can generate beams with OAM values of $=\pm 2$ or induce singular polarisation
patterns, e.g. vector vortex beams. OAM beams, as well as vector vortex beams,
have been exploited to encode classical and quantum information Nagali  _et
al._ (2009); D’Ambrosio  _et al._ (2012); Fickler  _et al._ (2012); Vallone
_et al._ (2014); Willner  _et al._ (2017); Sit  _et al._ (2018).
Spherulite-like structures are observed when evaporating solutions of ascorbic
acid dissolved in water (or other solvents, like ethanol). We confirmed our
prediction by performing interferometric experiments to ascertain the OAM
character of the transmitted beam and spatial Stokes polarimetry to
reconstruct the generated polarisation patterns during and after the
spherulite formation.

Figure 1: Images of ascorbic acid samples between cross polarizers,
illuminated by white light. A. Sample of ascorbic acid crystals grown at room
temperature, from a 10:1 ethanol:water solution. Many “maltese-cross”-like
patterns are visible, delimited by boundary regions. In the top left corner, a
singled out pattern is enlarged. B. Sample of ascorbic acid crystals grown at
$80^{\circ}$C, from a pure water solution. The “maltese-cross”-like patterns
exhibit a double-banded structure visible as concentric rings of alternating
colors. C. Sample grown in the same conditions as B. Point defects aligned
close together on a physical boundary cause a different pattern to arise in
the sample. The double-bandedness remains visible, without the cross pattern.
D. Sample of ascorbic acid crystals grown at $120^{\circ}$C, from a pure water
solution. A ”maltese-cross”-like pattern is visible, but appears disorganized
and coarse. The fibrils themselves are almost recognizable. The scale bar
corresponds to $\sim$0.5 mm.

##  II Results

Ascorbic acid ($C_{6}H_{8}O_{6}$), commonly referred to as vitamin C, forms
biaxial crystals with monoclinic spheroidal crystal units of $P2_{1}$
symmetry. Cox (1932); Hvoslef (1968). The unit cell and indicatrix ellipsoid
are shown in Fig. 2-A, and the molecular structure of ascorbic acid is shown
in Fig. 2-C. The molecules themselves are roughly planar and at room
temperature growth they are almost strictly oriented in the (010) plane (the
plane generated by the axes $a$ and $c$ in Fig. 2-A), the smallest refractive
index of the indicatrix is aligned in the $b$ direction.
In our experiment, commercially available ascorbic acid was spread over a
glass slide, dissolved in a solution which was made to evaporate and leave
free and dispersed ascorbic acid molecules in the sample. The choices of
solvent and temperature are important as they determine the thickness of the
sample (which affects the retardation), the rate of crystal formation, subtle
variations in crystal structure, and the amount of molecules present on the
slide (whether or not the effects will be visible). Water and ascorbic acid
have great affinity; however, due to its high surface tension water does not
spread out as well on the slide and can take up to a couple hours to dry, at
room temperature Uesaka and Kobayashi (2002). Therefore at first, we chose to
dissolve the ascorbic acid at room temperature in a solution of 10:1
ethanol:water Omar (2006). On the sample, the spherulite structures formed
around point defects from residual powder grains, either accidentally present
from undissolved powder grains in the solution or introduced manually after
evaporation. Uesaka and Kobayashi (2002); Ito  _et al._ (2003); Yamazaki  _et
al._ (2009). There is a large variety of patterns that can be observed under
different conditions of humidity and temperature. Yamazaki  _et al._ (2009);
Crist and Schultz (2016). In particular, in specific humidity conditions
(relative humidity $RH\approx 30-40\%$) and at room temperature Yamazaki  _et
al._ (2009), one can observe “maltese-cross”-like patterns when viewing the
sample between crossed polarizers and under white light illumination (see Fig.
1\- A). Under different conditions, for example at slightly higher
temperatures, the crystals can form with a twist, rotating around the
fibrillar axis at a certain angle, resulting in a periodic rotation of the
indicatrix ellipsoid along the radial direction Price (1959). This alternation
of adjacent regions of relative high and low birefringence forms “banded
maltese-cross”-like patterns, where, when viewed under cross-polarizers,
regularly spaced concentric rings of alternating brightness appear overtop of
the cross-like intensity pattern. In the case of ascorbic acid, “double-banded
maltese-cross” patterns form because of the biaxial nature of the rotating
ellipsoid [Fig. 1-B and C]. By choosing pure water as a solvent which we
heated at near boiling point, we observed the double banded patterns more
clearly. In more extreme conditions at higher temperatures, the spherulite
patterns become coarse and less organized (see Fig. 1-D).
In order to understand the action of ascorbic acid crystals on a light beam
around point defects we can adopt the model of a wave retarder with an
inhomogeneous optic axis orientation. This model works equally for biaxial
media once one fixes the propagation direction of the transmitted light. In
such a case there are still two (effective) refractive indexes associated with
specific orthogonal polarizations, which can be calculated knowing the
wavevector of the incident beam and the orientation of the index ellipsoid
(see e.g. Saleh and Teich (2019)). In the following we conventionally define
as “optic axis” the axis corresponding to the highest refractive index
experienced by the incident light beam. Assuming the sample sufficiently thin
and lying in the plane $z=0$, we have an optic axis angle $\theta$ which is a
function of the azimuthal angle and radius in cylindrical coordinates:
$\theta(\rho,\phi)$. The action of an anisotropic uniaxial medium can be
expressed by the following Jones matrix, written in circular polarization
basis:

$\displaystyle
U(\delta,\theta)=\cos\left(\frac{\delta}{2}\right)\begin{bmatrix}1&0\\\
0&1\end{bmatrix}+i\sin\left(\frac{\delta}{2}\right)\begin{bmatrix}0&e^{-2i\theta}\\\
e^{2i\theta}&0\end{bmatrix},$ (1)

where the parameter $\delta$ is the optical retardation of the medium,
$\delta=2\pi(n_{1}-n_{2})d/\lambda$, with $n_{1,2}$ the two effective
refractive indices, $d$ the sample thickness and $\lambda$ the wavelength.
Note that $\delta$ can be position dependent as well, and in particular, in
the case of banded “maltese-cross”-like patterns.

Figure 2: A. Schematic of the ascorbic acid crystal unit cell and indicatrix
ellipsoid. The three axes of the unit cell are: $a=17.299(8)$ Å, $b=6.353(3)$
Å, $c=6.411(3)$ Å, with $\omega=102^{\circ}11^{\prime}(08^{\prime})$ the angle
between $a$ and $c$. The three axes of the ellipsoid are:
$\alpha=1.46,\,\beta=1.6,\,\gamma=1.75$, with $\alpha$ parallel to $b$. B.
Twisting of indicatrix ellipsoid along the fibril growth directions. C. The
molecular structure of ascorbic acid: grey, red and white spheres represent
carbon, oxygen and hydrogen atoms respectively.

At half-wave retardation $\delta=\pi$, the effect of such a medium is to
invert the handedness of circular polarisation and introduce a phase factor
$\exp(\pm i2\theta(\rho,\phi))$, where the $+$ sign is in the case of input
left circular polarisation, and the $-$ sign for input right circular
polarisation. In particular, mediums characterized the optic axis angle
anisotropy

$\theta_{q}(\rho,\phi)=q\phi+\theta_{0},$ (2)

where $q$ is the topological charge of the molecular director field and
$\theta_{0}$ an offset angle, introduce azimuthal phase factors of $\exp(\pm
i2\theta_{0})\exp(\pm i2q\phi)$. Beams characterized by $\exp(i\ell\phi)$
factors are known for carrying a well-defined amount of OAM equal to
$N\,\hbar\ell$ (where $N$ is the average number of photons). Therefore, the
medium described by $U(\pi,\theta_{q})$ is converting the input spin angular
momentum (SAM) of each input photon into OAM, a process known as SAM-to-OAM
conversion (STOC), first investigated with $q$-plates Marrucci  _et al._
(2006b). Moreover, an input beam that is linearly polarized (or in a generic
elliptical polarisation state) will be transformed in a non-separable state of
polarisation and OAM which manifests in inhomogeneous polarisation patterns
Cardano  _et al._ (2012, 2013); D’Errico  _et al._ (2017).
For radial patterns, which we can expect from the radially directed fibrillar
structure of spherulites, one has $q=1$ and $\theta_{0}=0$. In general, since
the optic axis orientation is defined modulo $\pi$, $q$ can be either an
integer or a semi-integer number, $q=1/2$ being thus the smallest topological
charge that can characterize such defects. It is well known that, in nature,
topological defects typically manifest the smallest possible charge (since
this correspond to a smaller defect energy Chaikin and Lubensky (1995)).
However, in this case there are no constraints on the molecular inclination
along the $z$-axis, so semi-integer topological charges can still induce a
discontinuity in the pattern under a $2\pi$ rotation, as is typically observed
in umbilical defects Chandrasekhar (1992); Barboza  _et al._ (2015). To test
that our samples of ascorbic acid spherulites have a local optic axis around
defects of $\theta=\phi$, we observed both the STOC phenomenon and the
generation of inhomogeneous polarisation patterns.

Experimental characterization of the optic axis pattern. We first confirmed
that the ascorbic acid crystals orient radially by measuring the optic axis
pattern $\theta:=\theta(\rho,\phi)$. This can be done by shining the sample
with circular polarisation and projecting the transmitted beam on a linear
polarisation state. The resulting intensity distribution is given by

$I(\rho,\phi)\propto 1+\sin(\delta)\sin(2\theta(\rho,\phi)),$ (3)

where we assumed input right-circular polarisation and a projection on
horizontal polarisation. The measurement was performed with a 633 nm diode
laser source. At this wavelength our sample exhibits a retardation close to
$\delta=\pi/2$, ensuring a good visibility of the intensity modulation in the
measured pattern. As shown in Fig. 3-B, near point defects, the recorded
intensity distribution (averaged over different radial positions) is well
described by the relation $\theta=\phi$, associated with a radial distribution
of the molecular director around the defect. In proximity of boundaries which
separate regions containing point defects, we observe a jump in the optic axis
angle.

Figure 3: A. Schematics of the interferometry and polarimetry setups. The
reference beam’s polarization is rotated in the reference arm to ensure non-
orthogonal polarizations at the output of the BS. A lens (not shown) is also
added to the reference arm to induce a spherical wavefront onto the reference
gaussian beam. For polarization tomography, the reference beam is ignored and
the camera is placed between the PBS and the BS, as shown. B. Intensity
patterns characterizing the optic axis angle in two different regions of the
sample. Right-circular polarization is incident on the sample and projected
onto horizontal polarization, resulting in an intensity distribution related
to the optic axis angle by 3. On the right, the distribution around a defect.
It matches the simulated intensity for $\theta=\phi$, shown in the bottom
right corner. On the bottom left, the intensity distribution of the image
above, showing a region near a physical boundary. The intensity does not vary
in that region. The scale bar corresponds to 0.5mm. C. On top, interference
pattern generated with input right circular polarization, compared to the
theoretical pattern between a spherical gaussian and a beam carrying OAM +2.
On the bottom, interference pattern generated with input left circular
polarization, compared to the theoretical pattern between a spherical gaussian
and a beam carrying OAM -2. D. Density plot of the major axes of the point-by-
point polarization ellipses from a beam transmitted through a defect region of
the sample. The V-point is visible. A simulation of the expected density plot
is shown in the bottom right corner.

STOC in ascorbic acid samples. A simple way to ascertain the generation of OAM
in a laser beam is to look at its interference with a reference Gaussian beam
Marrucci  _et al._ (2006b). It is easy to show that, in the case of co-
propagating beams, the interference pattern consists of spiral fringes where
the number of arms is dictated by the OAM absolute value of the analyzed beam,
and the handedness by the OAM sign (see insets in Fig. 3-C for simulations).
We thus verified this result by inserting the samples in one arm of an
interferometer and recording the interference pattern for right handed and
left handed circular polarisation inputs (the experimental setup is shown in
Fig. 3-A). We used a laser beam at 525 nm, where the optical retardation of
the samples is typically close to $\pi$.
The resulting interference patterns are shown in Fig. 3-C for right and left
circularly polarisations, respectively, showing that the STOC phenomenon is
realized.

Polarimetry of vector modes. As the next step we illuminated the sample with
the same laser beam but with linear polarisation. The $q$-plate model predicts
that the resulting field, on the sample plane and at half-wave retardation, is
given by

$\mathbf{E}_{out}(\rho,\phi)\propto(e^{i2\phi}\mathbf{e}_{R}+e^{-i2\phi}\mathbf{e}_{L})e^{-\rho^{2}/w^{2}},$
(4)

where we assumed a Gaussian envelope field with waist $w$, and
$\mathbf{e}_{R,L}=\left(\mathbf{x}\pm i\mathbf{y}\right)/\sqrt{2}$, with
$\mathbf{x}$ and $\mathbf{y}$ the unit vectors of the $x$ and $y$ axes. Such a
field is linearly polarised everywhere but with an angle of the polarisation
ellipse given by $2\phi$ and a singular point on the beam axis (known as a
“V”-point) Cardano  _et al._ (2013); D’Errico  _et al._ (2017). We observed
the generation of these patterns by a point-by-point measurement of the Stokes
parameters, defined by: $s_{0}=\mathcal{I}_{total}$, $s_{1}=I_{H}-I_{V}$,
$s_{2}=I_{A}-I_{D}$, $s_{3}=I_{L}-I_{R}$, where $I_{H,V,A,D,L,R}$ are the
intensities of, respectively, horizontal, vertical, anti-diagonal, diagonal,
left circular and right circular polarisation, and $\mathcal{I}_{total}$
stands for total intensity. These intensities, and consequently the Stocks
parameters, can all be easily measured by a sequence containing a quarter-wave
plate, a half-wave plate and a polarizing beam splitter (Fig. 3-A). By
recording the different projection intensities on a camera, we reconstructed
the local polarisation ellipse. The resulting pattern, is qualitatively close
to the expected one. In particular the density plot of the polarisation
ellipse orientation (which is the phase of the field $s_{1}+is_{2}$) clearly
shows a singular pattern (Fig. 3-D). Due to the optical retardation not being
exactly $\pi$, the “V”-point is split in 4 points of circular polarisation
(“C”-points) each with topological charge 1/2 (see also ref. D’Errico  _et
al._ (2017)).
In Fig. 4, we show some frames of a recording of the Stokes field evolution
during the evaporation of the solvent, at room temperature. It can be seen how
inhomogeneous polarisation distributions start to develop around the defect
and to expand radially. Again, a V-point formation in the polarisation pattern
is clearly visible.

Figure 4: Frames showing the evolution of polarization distribution during the
growth of the sample at t=0min, 6.5min, 12.8min, 19.1min. The sample was grown
at room temperature, from a solution of 10:1 ethanol:water. a. Intensity of
the beam between cross polarizers at different times during crystal growth. b.
Plot of the major axis angle of the point by point polarization ellipses at
different times during crystal growth. The scale bar corresponds to $\sim$1
mm.

##  III Discussion and Conclusions

We have shown that ascorbic acid spherulites crystallize around defects, in
the absence of twisting, with a local optic axis pattern of $\theta=\phi$. We
have also shown experimentally the generation of OAM modes and vector beams.
Our findings show how ascorbic acid can be exploited as a material for
fabricating devices which exploit photonic spin-orbit phenomena.
The radial growth around the defect is likely the most energetically efficient
way to form the hydrogen bonds between atoms inside the crystal. It can be
interesting to explore the different crystallization structures when ascorbic
acid is grown on a polar substrate, instead of around a defect on a neutral
glass substrate.
Ascorbic acid is part of a large family of spherulites, each of which having
structures suggesting certain anisotropy of optic axis patterns around
defects. It would be interesting to explore the STOC phenomena in other kinds
of spherulites, which may have advantages over ascorbic acid. In particular,
recently the controlled generation of quartz spherulites Zhou  _et al._ (2021)
may offer an interesting way for creating vortex beams in the Tera Hertz
frequency domain Wang  _et al._ (2020) ($q$plates for this regime have been
recently realized using either 3D printing techniques Hernandez-Serrano  _et
al._ (2017) or quartz segmented waveplates Minasyan  _et al._ (2017)) It would
also shed light on the role of different crystal structures and symmetries in
the shaping of light.

##  IV Acknowledgement

This work was supported by the Canada Research Chairs (CRC) and Canada First
Research Excellence Fund (CFREF) Program, and Joint Centre for Extreme
Photonics (JCEP).

## References

  * Crist and Schultz (2016) Buckley Crist and Jerold M Schultz,  _Polymer spherulites: A critical review,_ Progress in Science of Polymers 56, 1 (2016).
  * Marrucci  _et al._ (2006a) Lorenzo Marrucci, C Manzo, and D Paparo,  _Pancharatnam-berry phase optical elements for wave front shaping in the visible domain: switchable helical mode generation,_ Applied Physics Letters 88, 221102 (2006a).
  * Marrucci  _et al._ (2006b) Lorenzo Marrucci, C Manzo, and D Paparo,  _Optical spin-to-orbital angular momentum conversion in inhomogeneous anisotropic media,_ Phys. Rev. Lett. 96, 163905 (2006b).
  * Rubano  _et al._ (2019) Andrea Rubano, Filippo Cardano, Bruno Piccirillo, and Lorenzo Marrucci,  _Q-plate technology: a progress review,_ JOSA B 36, D70 (2019).
  * Nagali  _et al._ (2009) Eleonora Nagali, Fabio Sciarrino, Francesco De Martini, Lorenzo Marrucci, Bruno Piccirillo, Ebrahim Karimi, and Enrico Santamato,  _Quantum information transfer from spin to orbital angular momentum of photons,_Phys. Rev. Lett. 103, 013601 (2009).
  * Cardano  _et al._ (2015) Filippo Cardano, Francesco Massa, Hammam Qassim, Ebrahim Karimi, S. Slussarenko, Domenico Paparo, Corrado de Lisio, Fabio Sciarrino, Enrico Santamato, Robert W. Boyd, and Lorenzo Marrucci,  _Quantum walks and wavepacket dynamics on a lattice with twisted photons,_Science Advances 1, e1500087 (2015), arXiv:arXiv:1407.5424v1.
  * Yan  _et al._ (2015) Lu Yan, P Gregg, Ebrahim Karimi, A Rubano, L Marrucci, R Boyd, and S Ramachandran,  _Q-plate enabled spectrally diverse orbital-angular-momentum conversion for stimulated emission depletion microscopy,_ Optica 2, 900 (2015).
  * Nivas  _et al._ (2015) Jijil JJ Nivas, Shutong He, Andrea Rubano, Antonio Vecchione, Domenico Paparo, Lorenzo Marrucci, Riccardo Bruzzese, and Salvatore Amoruso,  _Direct femtosecond laser surface structuring with optical vortex beams generated by a q-plate,_ Scientific reports 5, 1 (2015).
  * Karimi  _et al._ (2009) Ebrahim Karimi, Bruno Piccirillo, Eleonora Nagali, Lorenzo Marrucci, and Enrico Santamato,  _Efficient generation and sorting of orbital angular momentum eigenmodes of light by thermally tuned q-plates,_ Applied Physics Letters 94, 231124 (2009).
  * Piccirillo  _et al._ (2010) Bruno Piccirillo, Vincenzo D’Ambrosio, Sergei Slussarenko, Lorenzo Marrucci, and Enrico Santamato,  _Photon spin-to-orbital angular momentum conversion via an electrically tunable q-plate,_Applied Letters 97, 4085 (2010), arXiv:1010.4473.
  * Beresna  _et al._ (2011) Martynas Beresna, Mindaugas Gecevičius, Peter G Kazansky, and Titas Gertus,  _Radially polarized optical vortex converter created by femtosecond laser nanostructuring of glass,_ Applied Physics Letters 98, 201101 (2011).
  * Karimi  _et al._ (2014) Ebrahim Karimi, Sebastian A Schulz, Israel De Leon, Hammam Qassim, Jeremy Upham, and Robert W Boyd,  _Generating optical orbital angular momentum at visible wavelengths using a plasmonic metasurface,_ Light: Science & Applications 3, e167 (2014).
  * D’Ambrosio  _et al._ (2012) Vincenzo D’Ambrosio, Eleonora Nagali, Stephen P Walborn, Leandro Aolita, Sergei Slussarenko, Lorenzo Marrucci, and Fabio Sciarrino,  _Complete experimental toolbox for alignment-free quantum communication,_ Nat. Commun. 3, 961 (2012).
  * Fickler  _et al._ (2012) Robert Fickler, Radek Lapkiewicz, William N Plick, Mario Krenn, Christoph Schaeff, Sven Ramelow, and Anton Zeilinger,  _Quantum entanglement of high angular momenta,_ Science 338, 640 (2012).
  * Vallone  _et al._ (2014) Giuseppe Vallone, Vincenzo D’Ambrosio, Anna Sponselli, Sergei Slussarenko, Lorenzo Marrucci, Fabio Sciarrino, and Paolo Villoresi,  _Free-space quantum key distribution by rotation-invariant twisted photons,_ Phys. Rev. Lett. 113, 060503 (2014).
  * Willner  _et al._ (2017) Alan E Willner, Yongxiong Ren, Guodong Xie, Yan Yan, Long Li, Zhe Zhao, Jian Wang, Moshe Tur, Andreas F Molisch, and Solyman Ashrafi,  _Recent advances in high-capacity free-space optical and radio-frequency communications using angular momentum multiplexing,_ Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 375, 20150439 (2017).
  * Sit  _et al._ (2018) Alicia Sit, Robert Fickler, Fatimah Alsaiari, Frédéric Bouchard, Hugo Larocque, Patrick Gregg, Lu Yan, Robert W. Boyd, Siddharth Ramachandran, and Ebrahim Karimi,  _Quantum cryptography with structured photons through a vortex fiber,_Opt. Lett. 43, 4108 (2018).
  * Cox (1932) EG Cox,  _Crystalline structure of hexuronic acid,_ Nature 130, 205 (1932).
  * Hvoslef (1968) Jan Hvoslef,  _The crystal structure of l-ascorbic acid, ’vitamin c’. i. the x-ray analysis,_ Acta Crystallographica Section B: Structural Crystallography and Crystal Chemistry 24, 23 (1968).
  * Uesaka and Kobayashi (2002) H Uesaka and R Kobayashi,  _Pattern formation in the crystallization of ascorbic acid,_ Journal of crystal growth 237, 132 (2002).
  * Omar (2006) Waid Omar,  _Effect of solvent composition on crystallization process of ascorbic acid,_ Chemical Engineering & Technology: Industrial Chemistry-Plant Equipment-Process Engineering-Biotechnology 29, 119 (2006).
  * Ito  _et al._ (2003) Miho Ito, Machiko Izui, Yoshihiro Yamazaki, and Mitsugu Matsushita,  _Morphological diversity in crystal growth of l-ascorbic acid dissolved in methanol,_ Journal of the Physical Society of Japan 72, 1384 (2003).
  * Yamazaki  _et al._ (2009) Yoshihiro Yamazaki, Hiroki Yoshino, Machiko Izui, Yukiko Sato, and Mitsugu Matsushita,  _Humidity–temperature dependence of domain growth of ascorbic acid crystal,_ Journal of the Physical Society of Japan 78, 074001 (2009).
  * Price (1959) Fraser P. Price,  _The structure of high polymer sperulites,_ Journal of Polymer Science 37, 71 (1959).
  * Saleh and Teich (2019) Bahaa EA Saleh and Malvin Carl Teich,  _Fundamentals of photonics_ (john Wiley & sons, 2019).
  * Cardano  _et al._ (2012) Filippo Cardano, Ebrahim Karimi, Sergei Slusarenko, Lorenzo Marrucci, Corrado de Lisio, and Enrico Santamato,  _Polarisation pattern of vector vortex beams generated by q-plates with different topological charges,_Appl. Opt. 51, C1 (2012).
  * Cardano  _et al._ (2013) Filippo Cardano, Ebrahim Karimi, Lorenzo Marrucci, Corrado de Lisio, and Enrico Santamato,  _Generation and dynamics of optical beams with polarization singularities,_ Optics express 21, 8815 (2013).
  * D’Errico  _et al._ (2017) Alessio D’Errico, Maria Maffei, Bruno Piccirillo, Corrado de Lisio, Filippo Cardano, and Lorenzo Marrucci,  _Topological features of vector vortex beams perturbed with uniformly polarized light,_ Scientific reports 7, 1 (2017).
  * Chaikin and Lubensky (1995) P. M. Chaikin and T. C. Lubensky,  _Principles of Condensed Matter Physics_ (Cambridge University Press, 1995).
  * Chandrasekhar (1992) S. Chandrasekhar,  _Liquid Crystals_, 2nd ed. (Cambridge University Press, 1992).
  * Barboza  _et al._ (2015) R. Barboza, U. Bortolozzo, M. G. Clerc, S. Residori, and E. Vidal-Henriquez,  _Optical vortex induction via light-matter interaction in liquid-crystal media,_Adv. Opt. Photon. 7, 635 (2015).
  * Zhou  _et al._ (2021) Silang Zhou, Jordi Antoja-Lleonart, Pavan Nukala, Vaclav Ocelik, Nick R Lutjes, and Beatriz Noheda,  _Crystallization of geo2 thin films into alpha-quartz: From spherulites to single crystals,_ Acta Materialia 215, 117069 (2021).
  * Wang  _et al._ (2020) Honggeng Wang, Qiying Song, Yi Cai, Qinggang Lin, Xiaowei Lu, Huangcheng Shangguan, Yuexia Ai, and Shixiang Xu,  _Recent advances in generation of terahertz vortex beams and applications,_ Chinese Physics B 29, 097404 (2020).
  * Hernandez-Serrano  _et al._ (2017) AI Hernandez-Serrano, E Castro-Camus, and D Lopez-Mago,  _q-plate for the generation of terahertz cylindrical vector beams fabricated by 3d printing,_ Journal of Infrared, Millimeter, and Terahertz Waves 38, 938 (2017).
  * Minasyan  _et al._ (2017) Amalya Minasyan, Clément Trovato, Jérôme Degert, Eric Freysz, Etienne Brasselet, and Emmanuel Abraham,  _Geometric phase shaping of terahertz vortex beams,_ Optics letters 42, 41 (2017).