# Structured foraging of soil predators
unveils functional responses to bacterial defenses

Fernando W. Rossine  Gabriel Vercelli  Corina E. Tarnita  Department of
Ecology & Evolutionary Biology, Princeton University, Guyot Hall, Princeton,
NJ 08544, USA  Thomas Gregor  Joseph Henry Laboratory of Physics & Lewis-
Sigler Institute for Integrative Genomics, Princeton University, Princeton, NJ
08544, USA  Department of Stem Cell and Developmental Biology, CNRS UMR3738,
Institut Pasteur, 25 rue du docteur Roux, 75015 Paris, France

## Abstract

Predators and their foraging strategies often determine ecosystem structure
and function. Yet, the role of protozoan predators in microbial soil
ecosystems remains elusive despite the importance of these ecosystems to
global biogeochemical cycles. In particular, amoebae — the most abundant soil
protozoan predators of bacteria — remineralize soil nutrients and shape the
bacterial community. However, their foraging strategies and their role as
microbial ecosystem engineers remain unknown. Here we present a multi-scale
approach, connecting microscopic single-cell analysis and macroscopic whole
ecosystem dynamics, to expose a phylogenetically widespread foraging strategy,
in which an amoeba population spontaneously partitions between cells with
fast, polarized movement and cells with slow, unpolarized movement. Such
differentiated motion gives rise to efficient colony expansion and consumption
of the bacterial substrate. From these insights we construct a theoretical
model that predicts how disturbances to amoeba growth rate and movement
disrupt their predation efficiency. These disturbances correspond to distinct
classes of bacterial defenses, which allows us to experimentally validate our
predictions. All considered, our characterization of amoeba foraging
identifies amoeba mobility, and not amoeba growth, as the core determinant of
predation efficiency and a key target for bacterial defense systems.

## Introduction

Throughout biomes and across spatial scales — from the vast African savannas
to the minute rhizosphere microcosms — foraging strategies are key to how
predators shape species diversity and ecosystem function[estes2011trophic].
Studies of macroscopic ecosystems have revealed that different aspects of
predator foraging strategies — e.g., exploration behavior, dietary selectivity
— determine prey encounter rates and the effectiveness of predation deterrents
such as plant secondary metabolites[iason2006behavioral, foley2005plant] or
porcupine quills[katzner2015most]. Yet, microbial foraging strategies — in
particular those of soil bacterivore protozoans — remain understudied due to
unique challenges: the opacity of soil impedes in loco behavioral
observation[brethauer2020impacts], and metagenomic methods are hindered by the
substantial amplification biases in protozoan sequences[geisen2018soil] and
the low abundance of protozoan predators relative to their bacterial
prey[crowther2019global]. All told, we do not know what strategies soil
protozoan predators employ to explore their microscopic landscape, nor do we
understand the consequences of such foraging strategies to bacterial
population dynamics or to the effectiveness of bacterial defenses.

Thus, characterizing foraging strategies is a crucial missing link to a
predictive mechanistic understanding of the well-documented ecological
importance of soil protozoan predators, of which amoebae are the most
widespread and abundant[bates2013global, adl2005dynamics,
clarholm1981protozoan]. The amoebae are generalist predators that mold soil
communities by consuming and limiting bacterial
populations[clarholm1981protozoan, rosenberg2009soil, kreuzer2006grazing] all
the while altering bacterial taxonomic diversity[rosenberg2009soil,
kreuzer2006grazing]. Ultimately, amoebae consume and remineralize a
significant portion of the total soil microbial productivity[de1994modelling,
elliott1977soil], leading to increased nutrient cycling[bonkowski2004protozoa,
clarholm1985interactions]. But the consequences of microbial predation go far
beyond the soil itself, ranging from changes to biogeochemical cycles to the
emergence of pathogens: predator-driven increases of bio-available nitrogen
fertilize plants[bonkowski2004protozoa, kreuzer2006grazing,
clarholm1985interactions], remineralized carbon returns to the
atmosphere[bonkowski2004protozoa], and the interactions between amoebae and
fungi are implicated in the evolution of fungal
pathogenicity[albuquerque2019hidden].

Amoebae have a molecular machinery that enables prey search and handling
behaviors required for complex foraging strategies. Amoeba cells can perceive
and integrate various cues — mechanical forces[riviere2007signaling],
temperature[hong1983thermotaxis], and attractant/repellent chemical
gradients[bonner1969acrasin, phillips2012secreted] — leading to oriented cell
movement, often towards prey. Once amoebae encounter bacteria they deploy a
lineup of lectins that bind to specific bacteria surface carbohydrates,
leading to selective prey handling and consumption[farinholt2019microbiome].
When scaled up to thousands of amoebae, such multi-factorial behaviors may
lead to complex emergent features of the predator-prey interaction. Indeed,
studies have suggested that amoebae at higher densities consume bacteria more
efficiently, pointing towards benefits of collective prey-
handling[rubin2019cooperative]. The locomotive behavior of a collection of
cells may also lead to emergent ecological properties. For instance,
experimental and theoretical work on bacteria have shown that the interplay
between nutrient consumption and cell movement orientation creates chemical
gradients that determine how the bacterial population expands into uncolonized
regions[cremer2019chemotaxis, narla2021traveling]. Analogously, other studies
have suggested that the secretion of chemorepellents might play a role in the
expansion of soil amoebae colonies[phillips2012secreted, phillips2010roco].
Overall, describing how the behavior of individual amoeba cells results in
emergent foraging strategies is fundamental to understanding the efficiency
with which amoebae deplete their bacterial prey.

In this study, we identify the foraging strategies that soil amoebae adopt
when invading a spatially structured bacterial matrix and show that such
foraging strategies determine both the rates of bacteria consumption as well
as the effectiveness of different functional classes of bacterial anti-
predator defenses. We employ a cross-scale approach, characterizing both the
microscopic movement patterns of single amoebae within the bacterial matrix,
and the emergent macroscopic spatial patterns of predation. We first
investigate the foraging strategy of Dictyostelium discoideum, a well studied
bacterivore soil amoeba species with a complex life cycle comprised of a
unicellular foraging stage and a multi-cellular dispersal
stage[bonner2015cellular]. To assess the generality of the observed patterns,
we compare D. discoideum’s foraging strategy with that of other soil amoebae
species with distinct life cycles. In all investigated species we find cell
behavior differentiation and the occurrence of coordinated cell movement. We
build a mathematical model that integrates our microscopic and macroscopic
observations to make predictions about the macroscopic consequences of
disturbing individual cell behaviors. We test these predictions using chemical
deterrents that mimic two functional classes of bacterial defenses — predator
movement deterrents and predator growth rate deterrents. Altogether, we
present an experimental and theoretical framework that quantitatively bridges
individual amoeba behaviors and self-organized growth, predator-prey dynamics,
and microbial defenses.

## Results

### Macroscopic and microscopic features of bacterial matrix invasion

Figure 1: Macroscopic features of D. discoideum invasion: A, Schematic shows a
D. discoideum colony growing in a bacterial lawn illuminated by a blue LED
projector with EGFP-labelled E. coli emitting green light. Sectors show the
same Petri dish at $t_{0}=0$, $t_{1}=1.5$, and $t_{2}=3.5$ days after amoeba
inoculation. Dashed blue line denotes amoeba colony boundary and dashed red
line denotes area of bacterial exhaustion. False color denotes proportion of
remaining bacteria. B-E, Show data collected for 16 Petri dishes. B-C,
Aggregate profiles of remaining bacteria as a function of distance from the
colony boundary for times $t_{1}$ (B) and $t_{2}$ (C), each calculated from 40
profiles taken from each of the 16 dishes. Insets show detail of bacterial
accumulation along colony boundary. D, Time courses for the effective radii
for both the colony and the depleted bacteria zone. Inset shows estimated
values of $V_{x}$. E, Logarithm of consumed bacteria proportion as a function
of distance from colony boundary. Inset shows estimated values of
${L_{\frac{1}{2}}}$. For (D,E) transparent lines show time courses for
individual dishes and dashed lines show mean time course.

To elucidate the dynamics of how D. discoideum invades the bacterial matrix we
first set out to characterize the macroscopic features of such an invasion.
Notably, to ensure that our observations reflect the natural behavior of D.
discoideum cells, we use natural isolates, rather than the more commonly used
lab strains, known for their anomalous colony morphologies[fey2013one]. We
prepared Petri dishes with homogeneous lawns of eGFP-expressing E. coli cells
grown for two days until reaching the stationary phase. At the center of these
lawns we inoculated about 100 cells of D. discoideum . These amoeba cells
consumed bacteria, divided, and moved, thereby establishing outwardly
expanding colonies (Fig. 1-A).

Early on during the bacterial matrix invasion, a highly fluorescent ring
develops around the initial D. discoideum inoculation zone due to an
accumulation of bacteria. This luminous ring demarcates the boundary of the
feeding front beyond which amoeba cells are not present. Thus the luminous
ring is used to track the expansion of the nascent amoeba colony (Fig. 1-A).
We find that D. discoideum has colony expansion rates of $V_{x}=5.2\pm
0.3\,{\rm\mu m/min}$ on average. Notably, colonies attain such expansion rates
from the moment that the luminous ring becomes visible and before any
substantial bacterial consumption occurs (Fig. 1-D). This is surprising
because classical models for expanding populations of isotropically moving
cells predict that colony expansion rates should accelerate until the amoeba
feeding front stabilizes, which can only occur after bacteria are depleted at
the center of the dish[murray2007mathematical]. This discrepancy between
classical models and our observations suggests that anisotropic cell movement
might be an important factor in D. discoideum colony expansion.

To further characterize the feeding front, we traced radial fluorescence
intensity profiles, normal to the colony boundary (Fig. S1-A,B). Because
fluorescence is lost when bacteria are consumed by amoebae (Fig. S1-C,D),
these profile measurements quantify the proportion of consumed bacteria as a
function of the distance from the colony boundary. After a transient period
lasting up to 2 days, the bacteria consumption profile adopts a constant shape
consisting of three distinct zones: the D. discoideum colony boundary, which
is marked by a local accumulation of bacteria; a zone of exponential
consumption of bacteria; and a zone exhausted of bacteria (Fig. 1-B,C). To
characterize the exponential growth of the amoeba colony, we introduce the
halving length (${L_{\frac{1}{2}}}$) of the feeding front. It measures the
characteristic distance across the feeding front such that the proportion of
consumed bacteria falls by half. For D. discoideum we measured halving lengths
of ${L_{\frac{1}{2}}}=1.45\pm 0.08\,{\rm mm}$ (Fig. 1-E).

$V_{x}$ and ${L_{\frac{1}{2}}}$ provide a temporal and a spatial descriptor,
respectively, for the macroscopic characterization of amoeba feeding fronts:
$V_{x}$ reveals how the feeding front advances in time and ${L_{\frac{1}{2}}}$
reveals how the feeding front is shaped (Fig. S1-E,F). Importantly, $V_{x}$
and ${L_{\frac{1}{2}}}$ are not just descriptors of the feeding front, but are
also quantities that underlie ecological aspects of the predator-prey
interaction: as the bacteria consumption profile stabilizes, $V_{x}$ sets the
global bacteria consumption rate (i.e. ), while ${L_{\frac{1}{2}}}$ sets the
area over which amoeba cells compete for resources by defining the spatial
extent of the feeding front.

We next set out to derive the observed macroscopic features of the bacterial
matrix invasion from the underlying microscopic organization and dynamics of
amoeba cells at the feeding front. To this end, we cut out agar slabs around
the feeding front, and placed them under a confocal microscope. Infusion of
the bacterial lawn with fluorescence (i.e. fluorescein), allows us to track
amoeba cells as dark regions within the fluorescent background (Fig. 2-A,B).
While amoeba cells were present at all depths of the approximately
$100\,{\rm\mu m}$ thick bacterial matrix, most amoeba cells are located either
at its surface or at the interface between the agar gel and the bacterial
matrix. At the very boundary of the feeding front, marked by the
macroscopically observed luminous ring, amoeba cells accumulate — both at the
surface of the bacterial matrix but also at the interface with the agar gel —
forming a sharp band (Fig. 2-A). We consider two amoeba cell populations in
particular: cells around the boundary of the feeding front (referred to as
edge cells) and cells $2.5\,{\rm mm}$ behind the boundary (referred to as
inner cells). Tracking individual cells of these two populations uncovers two
strikingly distinct behavioral patterns: inner cells moved in no preferential
direction, whereas edge cells displayed markedly polarized trajectories,
moving preferentially away from the colony center (Fig 2 H,I). Moreover, we
find that with an average speed of $\bar{v}_{e}=3.00\pm 0.06\,{\rm\mu m/min}$
edge cells move over three times faster than inner cells with
$\bar{v}_{i}=0.81\pm 0.01\,{\rm\mu m/min}$ (Fig. 2-C).

Figure 2: Microscopic features of D. discoideum invasion: A, Region of the
feeding front imaged with a confocal microscope. Expanded region shows
tresholded confocal images showing D. discoideum cells at different depths of
the bacterial matrix, with cells shown in black at the surface of the
bacterial matrix and cells shown in light gray at the bottom. Each box
corresponds to one field of view along the feeding front, ranging from the
inner cells (dark blue) to the edge cells (light blue). B, Schematic of the
agar slab with the feeding front showing amoeba cells at different depths. C,
Median speed (box plots) and direction (ring plots) for individual amoeba
trajectories, for inner and edge cells, either embedded in the native feeding
front or isolated from the colony. Whiskers of the box plots denote the
$5^{th}$ and $95^{th}$ percentiles. D-G, Mean amoeba displacement (MD) and
position variance (S2???) are shown as a function of time for inner and edge
cells. Black points denote data, blue lines are bootstrapped linear
regressions (after a relaxation time) and white lines are the median
regressions. Insets show bootstrapped estimates of advection and diffusion
coefficients. H-I, Sample trajectories of inner (H) and edge (I) cells.

Polarized movement of individual amoeba cells suggests that they respond to
some environmental cue (e.g. a bacterially produced chemoattractant, a
chemorepellent produced by amoeba cells, or rheological changes to the
bacterial matrix caused directly or indirectly by the amoeba cells). Such cues
could in fact be responsible for the observed difference between edge- and
inner-cell speeds: in many chemotaxis systems both, cell polarization and
higher cell speed, arise via a common response to chemical
gradients[armitage1997bacterial, keller1971model]. However, it is unclear
whether differences in polarization and speed are indeed linked. To test in
our system whether cell speed differences persist in the absence of
environmental cues, we collect edge cells and inner cells separately, washed
them, and placed them in buffer-filled glass bottom dishes. As expected, as
cells settled and resumed their pseudopodial motion in this new environment
that is devoid of directional cues, no polarized cell movement could be
detected. However, edge cells remained faster than inner cells, indicating
that some of the behavioral differentiation between inner and edge cells can,
at least transiently, persist in the absence of environmental cues (Fig. 2-C).

The presence of individuals with high and low mobility in the same population
is a common pattern in ecology, showing up in species as disparate as sea
slugs[krug2009not], insects[roff1991wing], bacteria[armitage1997bacterial],
and even in the multicellular stage of the slime mold life
cycle[rossine2020eco]. This form of behavioral differentiation often reflects
the ubiquitous tradeoff between exploiting local resources and exploring the
landscape[mehlhorn2015unpacking]. Given the generality of this tradeoff, it is
natural to ask whether other species of soil amoebae would adopt similar
bacterial matrix invasion strategies. To investigate this question, we
selected another dictyostelid, Polysphondylium violaceum, and a more distantly
related amoeba species, Acanthamoeba castellanii, belonging to the Lobosa
group[cavalier2015multigene]. Despite being soil bacterivores, all three of
our species of interest are highly divergent and have very different life
cycles and cell characteristics. Surprisingly, we find that both P. violaceum
and A. castellanoo display the same general macroscopic and microscopic
invasion features that we uncovered in D. discoideum (Fig. S2). We observe in
all species the formation of a broad feeding front demarcated by a bacteria-
dense ring, a band of amoeba cells at the colony boundary, and behavioral
differentiation between slow-moving, unpolarized inner cells and fast-moving,
polarized edge cells.

### A microscopically-informed model for colony expansion

With both macroscopic and microscopic characterizations of the amoeba colony
expansion and bacterial matrix invasion at hand, we are now in a position to
ask whether our microscopic description of individual amoeba cell dynamics can
explain the macroscopic features of the invasion. In what follows we will use
our microscopic data to build a model that quantitatively predicts the
experimentally-measured macroscopic outcomes of the bacterial matrix invasion
and of the amoebal colony expansion.

We model radial amoeba density $\rho(x,t)$ as a continuous quantity, taking
into account cell division and movement, where $x$ is the position measured
relative to the amoeba colony boundary $b$. We assume amoebae divide at a
constant rate $r$ (given as a function of the mean doubling time DT) until
bacterial prey is locally exhausted, whereupon cell division halts. We
represent amoeba cell movement using a Fokker-Planck operator that consists of
an advection term to capture the polarized aspect of cell movement, and a
diffusive term to capture the dispersive aspect of cell movement (see SI). A
crucial difference between our model and classical taxis
models[keller1971model, cremer2019chemotaxis, narla2021traveling] our omission
of a mechanistic description of cell movement polarization. Instead, as
amoebae at different positions of the feeding front have different
characteristic movement behaviors (see above), we allow diffusion and
advection coefficients to be a function of the distance between the amoebae
and the colony boundary. We use our edge cell tracks to extract advection
($\nu_{e}$) and diffusion ($D_{e}$) values for cells around and beyond the
colony boundary, and our inner cell tracks to extract advection ($\nu_{i}$)
and diffusion ($D_{i}$) values for cells in the exponential consumption zone
of the feeding front (Fig. 2-D-I). Notably, given that inner cell movement is
unpolarized, the estimated value for inner cell advection is close to zero
(Fig. 2-D). To obtain advection and diffusion values for the transitional zone
(of length l) between the colony boundary and the exponential zone, we
interpolate the behavior of edge and inner cells. Following our observation
that inner and edge cells retain their differential behaviors even when
removed from the feeding front, we introduce a correlation time $\hat{t}$ for
which cells retain their movement pattern (see SI).

Assuming the existence of a stable invading wave solution, we find analytical
relations between the microscopic parameters of the model (i.e. $\nu_{e}$,
$\nu_{i}$, $D_{e}$, $D_{i}$, and DT) and the emergent macroscopic feeding
front features (i.e. $V_{x}$ and ${L_{\frac{1}{2}}}$):

$\displaystyle
V_{x}=\nu_{e}+2\cdot\sqrt{\frac{D_{e}\cdot\log(2)}{\textsf{DT}}}$ (1)
$\displaystyle{L_{\frac{1}{2}}}=\frac{V_{x}\cdot\textsf{DT}}{2}\cdot\left(1+\sqrt{1-\frac{D_{i}\cdot
4\log(2)}{V_{x}^{2}\cdot\textsf{DT}}}\right).$ (2)

Crucially, these analytical relations are independent of the details and
parameters (l, $\hat{t}$) of the behavioral interpolation across the
transitional zone. We can now employ these relations for an objective
criterion to evaluate the performance of our model: we simply require that
equations 1 and 2 hold for some combination of parameters within the ranges of
our measurements or literature-obtained values.

Figure 3: Modeling colony expansion: A-C, Results are shown for the complete
model with polarized edge cells (in blue) and for the model with unpolarized
edge cells (in red). A, Numerical realizations of the invading amoeba colony
profile produced by our model for the estimated ranges of $\nu_{e}$, $D_{e}$,
$D_{i}$, and a fit $K$. White lines show profiles for median estimated
parameters. Dashed black line shows data for bacterial consumption. B, Black
lined bell curve shows distribution of DT estimated from $D_{i}$, $V_{x}$, and
${L_{\frac{1}{2}}}$, and stripped area shows a literature range for DT. In
parallel $V_{x}$ is estimated from $\nu_{e}$, $D_{e}$, and DT. Black points
show measured values for $V_{x}$. C, Model predicted changes to $V_{x}$ both
in absolute terms (first row) or in percentage (second row) as we change cell
doubling time (first column), or cell speed (second column).

Using this criterion, we ask whether a single type of cell movement behavior
is able to account for the observed macroscopic features. We find that if
inner and edge cells all behave in the same way, either all as inner cells or
all as edge cells — which we implemented by setting $D_{e}=D_{i}$ and
$\nu_{e}=\nu_{i}$ — then the model cannot reconcile our macroscopic and
microscopic measurements (Fig. S3 A,B). However, when we integrate both of the
identified cell movement modes into our model, we find good agreement between
macroscopic and microscopic measurements (Fig. S3 C). This implies that the
macroscopic patterns of cell expansion cannot be attributed to any particular
cell movement mode. Rather, these patterns emerge from the interplay between
the distinct edge and inner cell behaviors.

Our model thus recapitulates a series of experimental observations. We find
that a narrow region of amoebae cells with polarized movement produce a broad
feeding front (Fig. 3-A), and our microscopic data accurately predicts the
expansion rate of this feeding front (Fig. 3-B). Moreover, by taking our cell
movement measurements and macroscopic colony measurements we employ equations
1 and 2 to estimate cell doubling times, which fall well within the range
found in the literature (Fig. 3-B). Surprisingly, despite the simplicity of
our assumptions about cell movement and cell division, the numerical solutions
of our model produce an accumulation of amoeba density at the boundary of the
colony (Fig. 3-A), which correspond to the sharp band of amoebae that we
experimentally identified.

To further investigate which factors of cell behavior lead to the formation of
the sharp amoeba band, we identified analytical conditions for the existence
of such a band in our model. The emergence of the invasion band can be
understood as a tug of war between edge cell movement and doubling times. Edge
cells with faster, more polarized motion — and therefore a higher $\nu_{e}$ —
favor the formation of the band by separating edge cells from the lagging
exponential zone, whereas faster doubling times counteract the separation of
the band by increasing cell density along the transitional zone (Fig. S4 A,B).
Parameters describing the transition between edge cell and inner cell behavior
also play a role in determining the formation of the invasion band. A narrower
transitional zone (smaller l) corresponded to a steeper transition between
edge and inner cell behavior, which favors the formation of the band.
Furthermore, a longer behavioral switching time $\hat{t}$ leads to increased
behavioral variance in the middle of the transitional zone, which also favors
the formation of the band (Fig. S4 C,D).

Finally, we used our model to investigate how disrupting cell behavior leads
to changes at the colony level. We consider three types of cell behavior
disruption. First, the effect of depolarizing the movement of edge cells,
which we find should lead to substantially decreased colony expansion rates
and narrower feeding fronts (Fig. 3-A,B). Second, we consider the effect of
changing doubling time, both by itself but also in conjunction with
depolarizing edge cell movement. Here we find that, by itself, lengthening of
doubling time should not substantially decrease colony expansion rates,
instead just broadening the feeding front (Fig. 3-A-C, Fig. S5). However, if
lengthening of doubling time is coupled to edge cell movement depolarization,
an appreciable decrease in colony expansion rate should occur. Third, we
consider the effect of changing cell speed (Fig. 3-A-C, Fig. S5), and find
that, regardless of edge cell polarization, colony expansion rates should be
directly proportional to cell speed (Fig. 3-C). Altogether these theoretical
results highlight the contribution of edge-cell polarization to the attained
colony expansion rate, as well as reveal that such movement polarization might
buffer the colony expansion rate from increased cell mortality and disrupted
cell division.

### Disrupting amoeba cell behavior

To probe our theoretical predictions of the colony-level consequences of
amoeba cell behavior disruption, we employed toxins that target distinct
aspects of amoeba behavior. First, we disrupt cell speed by using nystatin, a
toxin naturally produced by the soil-dwelling bacteria Streptomyces noursei.
Nystatin binds to sterols present in the cell membrane, a common eukaryote-
specific target for bacterial toxins, leading to altered cell membrane
properties such as increased permeability to ions and increased
rigidity[silva2006competitive]. Although nystatin causes no detectable
mortality to amoeba cells in our experimental conditions (Fig. S6), the speed
of individual edge cells is on average $25\pm 1\%$ slower in the presence of
nystatin compared to a control (Fig. 4-A). Consistent with the theoretical
predictions, the diminished cell speed in the experiments containing nystatin
results in a $25\pm 2\%$ slower expansion rate (Fig. 4-B) as well as in a
narrowing of the feeding front (Fig. 4-C,D). Second, we sought to
independently disrupt D. discoideum doubling time by slowing cell division. To
this end, we infused our bacterial lawns with fluorouracil, a synthetic uracil
analog that acts by inhibiting thymine synthesis to slow down DNA replication
and cell division[loomis1971sensitivity] but leaves the average edge cell
speed unaffected (Fig. 4-E). Consequently, consistent with our predictions, we
observe no reduction in the expansion rates of fluorouracil treated D.
discoideum colonies (Fig. 4-F), in conjunction with a broadening of the
feeding front (Fig. 4-G,H) due to the fluorouracil-induced lengthening of cell
doubling times. Overall, our toxin experiments corroborate the distinct
theoretically-predicted signatures of disrupting cell movement and doubling
time.

Figure 4: Toxin-disrupted D. discoideum behavior: A-D, Effects of nystatin
treatment. E-H, Effects of fluorouracil treatment. A,E, Changes to the median
speed of single amoeba trajectories for controls and toxin treatments.
Whiskers of the box plots denote the $5^{th}$ and $95^{th}$ percentiles. B,F,
Show the resulting differences to colony expansion rates. C,G, Representative
images of colonies growing under toxin treatments and controls. D,H, Aggregate
profiles of remaining bacteria as a function of distance from colony boundary,
both under toxin treatments and controls.

Toxins might also simultaneously affect both cell movement and doubling times.
Although nystatin did not reduce D. discoideum viability when cells were
embedded in the bacterial matrix (Fig. S6), high mortality is known to occur
for amoeba cells directly exposed to buffer with nystatin. Previous work has
shown that it is possible to select for D. discoideum variants that have lower
nystatin-induced mortality, and the resitance mechanism has been linked to
membrane sterol compositions with lower nystatin
binding[scandella1980genetic]. Since these resistant variants bind less
nystatin, we ask whether the motility effects of nystatin are also diminished
in these resistant cells, which should attenuate the reduction in colony
expansion rates caused by nystatin. We repeated our colony expansion
experiments with nystatin-resistant amoeba cells, which we selected following
previous work[scandella1980genetic]. Counterintuitively, we find that the
nascent colonies seeded from nystatin-survivor cells expand at a much slower
rate than those seeded from naive cells, even in the presence of nystatin
(Fig. 5-A,B). Moreover, our model predicts that such a reduction to the colony
expansion rate should be associated with the loss of the distinct invasion
band and a strong steepening of the feeding front (Fig. S4 A, Fig. 3 A), which
indeed we verify experimentally (Fig. 5-C,D).

The slowly expanding, nystatin-survivor phenotype cells are potentially at a
selective disadvantage when growing alongside naive cells. Accordingly, in the
absence of nystatin, we expect a quick loss of resistance. Indeed, a few days
into the expansion of colonies seeded with nystatin-survivor cells, small
sectors displaying the naive expansion pattern emerge from points around the
colony boundary. These naive-like sectors recover the formation of the broad
feeding front with the sharp invasion band, and eventually take over the
entire colony boundary (Fig. 5-D,E). Subsequently repeating the colony
expansion experiments with cells sampled from the naive-like sectors leads to
expansion patterns that are indistinguishable from those of colonies seeded
with truly naive cells (Fig. 5-A), suggesting a phenotypic reversal to the
naive state.

Figure 5: Properties of nystatin-tolerant colonies: A, Time courses of the
effective radius of the bacterial exhaustion zone for D. discoideum colonies
started from cells that were plated immediately after being selected for
nystatin-survival (early S cells, red), nystatin-naive (N) cells grown without
nystatin (solid blue lines), N cells grown with nystatin (solid gray lines),
and cells selected for nystatin-survival that were allowed to grow in
nystatin-infused Petri dishes for three days (S cells) grown with (dashed gray
lines) or without nystatin (dashed blue lines). B, Colony expansion rates for
first two days of growth. C-E, representative images of early colony growth
for N cells (C), early S cells (D) and the emergence of naive-like cells in a
colony founded by early S cells (E, same colony as D, but a day later).

Altogether, our results indicate that not only does targeting distinct aspects
of cell function lead to distinct responses at the colony level, but that
tradeoffs can emerge between different aspects of toxin tolerance. Phenotypes
that seemingly provide a viability advantage against a toxin might actually
underperform when other aspects of cell function, such as motility, are
considered.

### Discussion

When foraging within the bacterial matrix—a three dimensional environment with
many structural components—the soil-dwelling amoeba D. discoideum adopts a
spatially-structured foraging strategy with faster, polarized edge cells and
slower, unpolarized inner cells. This movement strategy gives rise to fast-
expanding colonies that develop a sharp ring of explorer cells. Several
distantly related soil-dwelling amoeba species all adopt the same general
structured movement mode, leading to colonies with ring-like expansion.
Curiously, similar ring-like expansion has been identified in phylogenetically
and structurally distinct systems, such as bacteria growing in culture
media[cremer2019chemotaxis]. The pervasiveness of this ring expansion pattern
— which occurs not only across distinct genera of amoebae but also across
different domains of life (Eukaryota and Bacteria), regardless of whether
cells move using pseudopods or beating flagella, over substrates as distinct
as culture media or a bacterial matrix — might be a product of both the
ubiquitous tradeoff between exploiting local resources and exploring the
environment[mehlhorn2015unpacking] and the energetic efficiency of the
underlying spatially-differentiated movement modes.

High motility edge cells ensure that the colony will attain a high $V_{x}$,
while allowing the low motility inner cells to conserve energy without
reducing the colony $V_{x}$. In the case of amoeba cells, the tradeoff between
cell motility and energy efficiency is even more acute than what would be
expected solely due the direct energetic costs of locomotion. Because the same
cytoskeletal machinery is involved both in pseudopod-based locomotion and in
endocytosis, faster moving cells have a slower bacteria intake and accrue less
biomass[chubb2000dictyostelium]. An interesting possibility is that the
behavior differences that we observed between inner and edge cells might arise
simply by changing the regulation of endocytosis and pseudopod formation. RasS
is a protein that is known to be involved in endocytosis regulation, and
indeed, D. discoideum RasS knockouts were shown to have higher rate of
pseudopod formation relative to endocytosis rate[chubb2000dictyostelium].
These mutant lineages display cell speed and polarized morphology very similar
to the ones we identified in our wild-type isolated edge cells, suggesting
that RasS might be involved in the behavioral differentiation between edge and
inner cells. Overall, this observed amoeba cell behavioral differentiation may
provide the means to attain a high $V_{x}$ while limiting the compromise to
bacteria consumption.

The emergent colony expansion rate $V_{x}$ could play an important ecological
role in the competition among amoebae for bacterial prey. Amoebae inhabit a
dynamic soil environment, in which foraging can only occur during short
metabolically active periods following rain or nutrient
pulses[clarholm1981protozoan, adl2005dynamics]. To persist, amoebae rely on
drought resistant spores that survive and disperse to new bacteria patches
during quiescent periods. Amoeba populations compete to produce and disperse
these spores. $V_{x}$ is likely to affect spore production in two ways:
modulating resource consumption and interference with competitors. Firstly, a
higher $V_{x}$ could allow faster-expanding amoebae to consume more bacteria
and thus produce more spores (Fig. S7-A,B). Secondly, if amoebae co-occur in
the same bacteria patch, the ones with a higher $V_{x}$ may deprive the lower
$V_{x}$ amoebae of these bacterial resources: faster-expanding amoebae
physically displace slower-expanding amoebae at the colony boundary
[buttery2012structured]. This prevents the slower amoebae from accessing
bacteria and producing spores (Fig. S7-A,B). Fast expansion (high $V_{x}$)
thus both increases one’s own dispersal capacity and may also reduce the
dispersal capacity of others. The relative contributions of these two
competitive modes to overall competition depend on ecological factors,
including how often competing types co-occur, the durations of metabolically
active and quiescent periods, and the size of bacteria patches. Nevertheless,
a higher $V_{x}$ should always produce a dispersal advantage: higher relative
bacteria consumption, more spore production, and consequently higher rates of
patch colonization (Fig. S7-A,B).

From the perspective of the bacterial prey, the ability to resist predator
invasion throughout the metabolically active soil periods might result in more
bacterial cells that can persist throughout the quiescent periods (spores or
otherwise), which would provide a population advantage when conditions improve
at the start of a new growth cycle (Fig. S7-C). Accordingly, bacteria deploy a
series of defenses against invading amoebae: structural elements of the
bacterial matrix may hamper the advance of invading
amoebae[huws2005protozoan], changes to bacteria cell shape impede
endocytosis[jurgens2002predation], and many toxic secondary metabolites are
synthesized[raaijmakers2012diversity, mazzola2009protozoan, klapper2018role,
klapper2016bacterial]. These defenses might act by sterilizing the amoebae and
preventing the invasion altogether, by reducing $V_{x}$ leading to a decrease
in bacteria consumption, or both. However, because the benefit of
sterilization without reduction of $V_{x}$ is substantially conditioned on
total predator elimination, such defense strategies are fragile to predator
adaptation or evolutionary change. Point in case, toxin tolerance varies
greatly amongst different protozoan predators in the same
ecosystem[brock2016sentinel, klapper2018role], evolution of resistance often
occurs in the face of previously sterilizing toxin
doses[scandella1980genetic], and toxic secondary metabolites are frequently
found in the soil in concentrations that are substantially lower than the
experimentally characterized minimum inhibitory
concentration[raaijmakers2012diversity, jousset2008secondary]. This suggests
that reduction of $V_{x}$ is an important component of bacterial defenses, and
our results imply that bacterial defenses that reduce amoeba doubling time —
either by increasing amoeba death rates or decreasing amoeba division rates —
lead to a limited reduction of $V_{x}$, whereas defenses that target cell
motility but that are not necessarily lethal to the predators — as we have
shown to potentially be the case for nystatin — achieve a more substantial
reduction of $V_{x}$.

If, from an ecological perspective, foraging strategies play an important role
in microbial predator-prey dynamics, from an evolutionary perspective these
same foraging strategies might have been pivotal to major evolutionary
transitions such as the emergence of multicellularity. In particular,
dictyostelids have a complex life cycle in which single cells forage until
food exhaustion at which point cells aggregate and develop into multicellular
fruiting bodies. When first deciphering the mechanisms that allow for cellular
aggregation in D.discoideum, John Bonner proposed that the cellular machinery
responsible for tracking food sources might have been repurposed into the
multicellularity pathway[bonner1969acrasin, bonner2015cellular]. Although it
was later established that dictyostelid signalling pathways were not directly
derived from prey chemosensing pathways[shimomura1982chemical], there is still
merit to the general idea that foraging strategies might have provided
preadaptations for the emergence of multicellularity. Here we have shown that
two important features of multicellular development, namely the possibility of
cell behavior differentiation and the occurrence of coordinated cell movement,
are present in the unicellular foraging phase of the dictyostelid life cycle.
Moreover, these features are also present in feeding fronts of Acanthamoebae,
a group that diverged from dictyostelids much earlier than the evolution of
the multicellular life cycle. This suggests that some pieces necessary for the
coordination of multicellular development are ancestral features of
amoebozoans, and helps shed light on the recurrent and independent appearances
of multicellular life forms across Amorphea, including metazoans, fungi, and
slime molds.

## Acknowledgements

We thank J. Bos for helpful comments and suggestions at various stages of the
work. FWR and CET acknowledge support from NSF RoL: FELS: EAGER-1838331. This
work was supported in part by the U.S. National Science Foundation, through
the Center for the Physics of Biological Function (PHY-1734030), and by
National Institutes of Health Grant R01GM097275.

Supplementary Information for
Structured foraging of soil predators
unveils functional responses to bacterial defenses

###### Contents

  1. A Supplemental Figures
  2. B Methods
    1. B.1 Cultures and strains
    2. B.2 Toxins
    3. B.3 Macroscopic imaging
    4. B.4 Analysis of colony boundaries
    5. B.5 Microscopic imaging
    6. B.6 Analysis of cell trajectories
    7. B.7 Cell viability
    8. B.8 A model for colony spread
    9. B.9 Mathematical relations between model quantities

##  Appendix A Supplemental Figures

Figure S1: Measuring the amoeba feeding front: A, white lines show location of
the profiles that were aggregated to construct the average bacteria
consumption profiles. B Aggregate bacteria consumption profiles for each of 16
independent Petri dishes. C Outer zone beyond the D. discoideum colony
boundary and depleted zone in which bacteria have been completely consumed. D
Normalized light intensity for pixels in the outer zone shows that there is no
artifact of irregular illumination or photobleaching. The low light intensity
across the depleted zone shows that any D. discoideum autofluorescence is
negligible, and that dugestion of bacteria effectively eliminates EGFP
fluorescence. E Closeup showing the colony boundary, the feeding front, and
the macroscopic quantities $V_{x}$ and ${L_{\frac{1}{2}}}$. F Effect of
different potential values of ${L_{\frac{1}{2}}}$ on the broadness of the
feeding front.

Figure S2: Macroscopic and macroscopic features of A.castellanii and P.
violaceum: A,B, Images of feeding fronts for P. violaceum (A) and A.
castellanii (B). False color denotes proportion of remaining bacteria. C,D,
Aggregate profiles of remaining bacteria as a function of distance from the
colony edge for P. violaceum (C) and A. castellanii (D), each calculated from
40 profiles taken from each of $15$ dishes. Insets show detail of bacterial
accumulation along colony edge. E,F, confocal microscopy views of the colony
edge for P. violaceum (E) and A. castellanii (F). G, Median speed and
direction (ring plots) for individual amoeba trajectories, for inner and edge
P. violaceum and A. castellanii cells. Whiskers of the box plots denote the
$5^{th}$ and $95^{th}$ percentiles.

Figure S3: Model agreement with data: A,B,Black lines show the theoretically
predicted relations between macroscopic and microscopic parameters if all
cells moved like inner cells (A) or if they all moved like edge cells (B.
Copper regions show the experimentally measured macroscopic quantities
($V_{x},{L_{\frac{1}{2}}}$) on the x axis, and the theoretically required
matching quantities on the y axis. Teal regions show experimentally measured
microscopic quantities on the y axis and the the theoretically required
macroscopic quantities on the x axis. The mismatch between copper and teal
regions means that under a homogeneous movement model the macroscopic and
microscopic data are inconsistent. C Incorporating differential movement rules
for inner and edge cells, heatmaps show theoretically predicted doubling times
as a function of $V_{x}$, ${L_{\frac{1}{2}}}$, and $D_{i}$. Star shows modal
experimental values, and region between dashed lines is the parameter region
that leads to doubling times between 2.7 and 5.5 hours (the literature range)

Figure S4: Conditions for the emergence of invasion band: A,C, region of the
parameter space such that a distinct band forms at the edge of the expansion
front for fixed $\hat{t}$ and $l$ and varying $\nu_{e}$ and DT (A) and for
fixed $\nu_{e}$ and DT and varying $\hat{t}$ and $l$ (C). B,D, numerical
realizations of the amoeba density around the colony edge for parameter
combinations color coded in A and C (B and D respectively).

Figure S5: Predicted effects of changing cell speed and doubling time:
Heatmaps show how jointly changing cell speed and doubling time leads to
theoretical changes in halving length and colony expansion rate. Results are
shown both incorporating polarized cell movement and in the absence of
polarized cell movement.

Figure S6: D. discoideum cell viability under nystatin treatment: A, honeycomb
lattice cut, nystatin infused Petri dish with D. discoideum colonies
independently growing from single D. discoideum cells in each hexagonal
lattice unit. B, number of lattice units colonized divided by expected number
of lattice units that received at least one viable D. discoideum cell for
control Petri dishes and nystatin-infused Petri dishes. Each point is the
aggregate count of viable D. discoideum colonies in $5$ Petri dishes with $48$
lattice units each. C likelihood profile of the effect of the addition of
nystatin on D. discoideum cell viability. Note that the maximum likelihood
value is very close to $1$ and that the likelihood ratio interval (shown in
blue) contains $1$, denoting no detectable effect of nystatin on D. discoideum
cell viability.

Figure S7: Ecological consequences of varying $V_{x}$: A, Schematic of a
bacteria patch being consumed by high $V_{x}$ and low $V_{x}$ amoebae. Upper
arrow points to relative amount of spores of each amoeba type produced after
the growth period, which is proportional to the bacteria area consumed by each
amoeba type. B, Schematic for relative amoeba spore production after long
growth cycles (shown in three time steps), in which amoebae consume all
bacteria, and after short growth cycles, in which some bacteria cells persist
after the growth period ends. Upper row shows relative spore production if
different amoebae types occur across separate patches (No interference
competition). Lower row shows relative spore production of amoebae co-
occurring in the same resource patch (interference competition is possible).
Note that, if different amoebae co-occur, long growth cycles maximize relative
spore differences, whereas, if different amoebae do not co-occur, long growth
cycles negate relative spore differences.

##  Appendix B Methods

###  B.1 Cultures and strains

Experiments used D. discoideum strain NC105.1 sourced from dictybase, A.
castellanii strain ATCC 30011 from ATCC, and a natural, locally collected P.
violaceum isolate. All amoebae were grown in E. coli lawns in adapted SM $2\%$
agar dishes with only $2.5g/L$ of dextrose. The pH of the agar dishes was
measured daily throughout experiments to ensure that it was always between 6
and 7. For imaging, agar had $2\%$ in mass of Mars black pigment (Gamblin) to
decrease gel fluorescence, added in prior to autoclaving. The E. coli carried
the EGFP expressing plasmid pEB1-mEGFP (a gift from Philippe Cluzel via
Addgene plasmid 103979 ; http://n2t.net/addgene:103979 ; RRID:Addgene_103979).
Homogeneous E. coli lawns were prepared by creating a mist of bacterial
culture using an air humidifier and letting the mist settle on the Petri
dishes.

###  B.2 Toxins

Nystatin in DMSO was added to culture media after autoclaving when the
temperature fell to 50 degrees celsius to avoid nystatin decomposition.
Nystatin concentration was 100 $\mu g/mL$ in treatment samples and DMSO
concentration was 2.5 $mL/L$ for both treatment and controls. For the
fluorouracil treatments, E. coli lawns were grown and then a fluorouracil
aqueous mist was let to settle over the agar gel for a final fluorouracil
concentration of 100 $\mu g/mL$. It is essential to add fluorouracil after the
bacteria lawn has grown, otherwise the drug will inhibit bacterial growth.

###  B.3 Macroscopic imaging

Our fluorescent illumination apparatus consisted of an AAXA P4-X Pico led
projector, supplemented with Roscolux CalColor 90 blue cinegel to reduce light
bleedthrough. Images were acquired by a DSLR EOS Cannon camera equipped with a
Canon RF 24-105mm f/4-7.1 IS STM lens enhanced with a nightsea 600nm emission
filter. At each experimental run, up to 30 Petri dishes were sequentially set
on an automated machine (built using a Welter’s rectangular lazy Susan as a
base) that took photos of each dish at 30 min intervals for a week.

###  B.4 Analysis of colony boundaries

Colony boundaries were determined by leveraging the existence of a highly
fluorescent ring around the colony. Microscopic inspection was used to
experimentally verify that the ring overlapped with the colony boundary.
During image analysis, the most fluorescent closed circumference on the dish
was found and its inner area taken to be the amoeba colony. The effective
radius of the colony was defined as the radius of a circle with the same area
as the colony.

###  B.5 Microscopic imaging

For microscopic imaging, commercial fluorescein-based water tracing dye
(EcoClean solutions) was added to media before autoclaving (10 m$L/L$). After
Amoeba colonies had grown for about three days, a slice of the gel around the
feeding front was removed and placed in a humidified chamber for imaging under
a Zeiss LSM microscope. Images were taken for up to 45 minutes after the
collection of the slice.

###  B.6 Analysis of cell trajectories

Individual cell trajectories were extracted after image processing using the
track mate plugin from Fiji. Confidence intervals for Diffusion and advection
coefficients were calculated by bootstrapping: trajectories were resampled
with substitution to re-calculate the coefficients.

###  B.7 Cell viability

To test the effects of nystatin on cell viability we laser cut an hexagonal
lattice on the surface of a nutrient agar gel, after bacterial misting, but
before lawn growth. After bacterial lawn growth, we added 5 $\mu L$ of D.
Discoideum cell suspension (at a concentration of approximately one cell per
droplet) to each lawn hexagon. Three days later we photographed the Petri
dishes and counted the number of viable colonies in treatment versus control.
For parameter estimation, colony counts were modelled as a poisson
distribution, and cell viability was modelled as a multiplier to the poisson
distribution for experiments with nystatin.

###  B.8 A model for colony spread

We model amoebae density $\rho$ as a continuous quantity that changes with
amoeba growth and movement. We assume that amoeba behavior is a function of
cell position relative to the amoeba colony boundary ($b$). Because of the
continuity of our density description, we define $b$ as being the last point
in space that has an amoeba density larger than some arbitrary fraction $f$ of
a carrying capacity $K$:

$\displaystyle b(t)=sup\\{x\in\mathbb{R}\mid\rho(x,t)>f\cdot K\\}$ (3)

Putting together growth and movement, we have that

$\displaystyle\frac{\partial}{\partial
t}\rho(x,t)=F_{x}(x,t,\rho)+r\cdot\rho(x,t)\cdot(1-C(K,\rho))$ (4)

Where $F_{x}(x,t,\rho)$ is a Fokker-Planck operator that incorporates the
randomness of amoeba movement into an diffusion coefficient $D_{x}(x-b)$ and
the biases of amoeba movement into an advection coefficient $\nu_{x}(x-b)$:

$\displaystyle F_{x}(x,t,\rho)=-\frac{\partial}{\partial
x}\left[\nu_{x}(x-b)\cdot\rho(x,t)\right]+\frac{\partial^{2}}{\partial
x^{2}}\left[D_{x}(x-b)\cdot\rho(x,t)\right]$ (5)

We assume that amoebae divide at a constant rate $r$ as long as they have
bacteria to prey upon:

$\displaystyle C(K,\rho)=\begin{cases}0&\text{if $\rho(x,t)<K$}\\\ 1&\text{if
$\rho(x,t))\geq K$}\end{cases}$ (6)

We use experimentally measured values of $\nu_{x}$ and $D_{x}$ for edge cells
($\nu_{e},D_{e}$) and for inner cells ($\nu_{i},D_{i}$). Following our
observation that inner and edge cells retain their differential behaviors even
when removed from the feeding front, we allow cells to retain their movement
patterns for a certain correlation time $\hat{t}$, after which they reevaluate
their behaviors as a function of their position relative to the boundary
($\bar{x}=x-b$). For simplicity, we say that the probability of a cell
behaving as an edge cell decreases linearly from 1 at the colony boundary to 0
at a point inside of the feeding front at a distance $l$ from the boundary
($p_{e}=1+\bar{x}/l$). This leads to advection and diffusion coefficients
given by

$\displaystyle\nu_{x}(\bar{x})=\begin{cases}\nu_{i}&\text{if $\bar{x}<-l$}\\\
\nu_{e}&\text{if $\bar{x}>0$}\\\ \nu_{i}\cdot(1-p_{e})+\nu_{e}\cdot
p_{e}&\text{if $-l\leq\bar{x}\leq 0$}\end{cases}$ (7) $\displaystyle
D_{x}(\bar{x})=\begin{cases}D_{i}&\text{if $\bar{x}<-l$}\\\ D_{e}&\text{if
$\bar{x}>0$}\\\ D_{i}\cdot(1-p_{e})+D_{e}\cdot
p_{e}+\hat{t}\cdot\frac{p_{e}\cdot(1-p_{e})\cdot(v_{e}-v_{i})^{2}}{2}&\text{if
$-l\leq\bar{x}\leq 0$}\end{cases}$ (8)

Note that the advection coefficient simply linearly interpolates the inner and
edge values. However, the diffusion coefficient may have an intermediate
maximum. This occurs because within the transitional zone there is a mix of
cells behaving either like inner or edge cells, in other words, movement is
more unpredictable in the transitional zone. However, if $\hat{t}$ is 0, then
cells have no memory and intermediary diffusion is just the linear
interpolation.

Importantly, analogously to the derivation of Fisher’s equation, we do not
derive our equations directly from the integrated stochastic cellular movement
and division processes. Instead, we simply add together growth and movement
terms that were independently derived. This procedure is known to lead to
pathological behaviors for extreme parameter values — in the case of Fisher’s
equation, for instance, if growth rates are high enough, one might obtain
travelling wave solutions that are faster than the maximum individual speed,
which is physically impossible. However, such pathological behaviors only
occur for parameter regions far from the experimentally measured quantities,
and indeed, individual based simulations indicate that the equations are
adequately representing the underlying stochastic process for the parameter
range of interest.

###  B.9 Mathematical relations between model quantities

The standard Fisher’s equation — which differs from our equation in having
constant diffusion and advection terms — has a unique stable travelling wave
solution that is asymptotically achieved from initial conditions with compact
support (that is, there are infinite possible travelling wave solutions, but
only one can be attained if at an initial moment all population density is
constrained to a finite region). We do not prove that this property holds for
our modified equation, but all of our numerical investigations suggest that
this is the case. In this section we assume that this is true and use this
property to calculate the $V_{x}$ and ${L_{\frac{1}{2}}}$ associated to this
asymptotic travelling wave solution.

Since the dynamics around and beyond the colony boundary follow Fisher’s
equation with constant diffusion and advection coefficients, and since it is
the front of the wave that determines wave speed in such cases, following the
classical derivation presented by Murray[murray2007mathematical] the colony
expansion speed $V_{x}$ must be given by equation 1.

To find ${L_{\frac{1}{2}}}$ we separately solve for the back of the wave.
Assuming $K$ is large enough, that there is indeed a wave solution with speed
$V_{x}$, and considering that $\nu_{i}=0$, then, for $\bar{x}<-l$ equation 4
becomes

$\displaystyle
D_{i}\frac{\partial^{2}}{\partial\bar{x}^{2}}\rho(\bar{x})+V_{x}\frac{\partial}{\partial\bar{x}}\rho(\bar{x})+r\cdot\rho(\bar{x})=0$
(9)

Which has an exponential solution. Since ${L_{\frac{1}{2}}}$ is just $\log 2$
divided by the exponent of $\rho$, it can be easily solved leading to equation
2.

Finally, to find the conditions for the emergence of a separate invasion band,
we note that all we require is that the left hand derivative of
$\rho(\bar{x})$ calculated at 0 be positive. In other words, separation occurs
if the amoeba density increases as we approach the colony boundary. A few
calculations substituting equations 8 and 7 into 4 lead to the condition that

$\displaystyle r<\hat{t}\cdot\frac{\nu_{e}}{l}+\frac{\nu_{e}}{l}$ (10)


# Magnetoresistance and Kondo effect in the nodal-line semimetal VAs2

Shuijin Chen  Zhefeng Lou  Qin Chen  Binjie Xu  Chunxiang Wu  Department of
Physics, Zhejiang University, Hangzhou $310027$, China  Jianhua Du  Department
of Applied Physics, China Jiliang University, Hangzhou $310018$, China  Jinhu
Yang  Department of Physics, Hangzhou Normal University, Hangzhou $310036$,
China  Hangdong Wang  Department of Physics, Hangzhou Normal University,
Hangzhou $310036$, China  Minghu Fang  Corresponding author<EMAIL_ADDRESS>Department of Physics, Zhejiang University, Hangzhou $310027$, China
Collaborative Innovation Center of Advanced Microstructures, Nanjing
University, Nanjing $210093$, China

###### Abstract

We performed calculations of the electronic band structure and the Fermi
surface as well as measured the longitudinal resistivity $\rho_{xx}(T,H)$,
Hall resistivity $\rho_{xy}(T,H)$, and magnetic susceptibility as a function
of temperature and various magnetic fields for VAs2 with a monoclinic crystal
structure. The band structure calculations show that VAs2 is a nodal-line
semimetal when spin-orbit coupling is ignored. The emergence of a minimum at
around 11 K in $\rho_{xx}(T)$ measured at $H$ = 0 demonstrates that an
additional magnetic impurity (V4+, $S$ = 1/2) occurs in VAs2 single crystals,
evidenced by both the fitting of $\rho_{xx}(T)$ data and the susceptibility
measurements. It was found that a large positive magnetoresistance (MR)
reaching 649% at 10 K and 9 T, its nearly quadratic field dependence, and a
field-induced up-turn behavior of $\rho_{xx}(T)$ emerge also in VAs2, although
MR is not so large due to the existence of additional scattering compared with
other topological nontrival/trival semimetals. The observed properties are
attributed to a perfect charge-carrier compensation, which is evidenced by
both calculations relying on the Fermi surface and the Hall resistivity
measurements. These results indicate that the compounds containing V
($3d^{3}4s^{2}$) element as a platform for studying the influence of magnetic
impurities to the topological properties.

###### pacs:

Since the discovery of quantum Hall effect in the two-dimensional electron gas
in the 1980s Laughlin (1983), it has been recognized that the topology of band
structure in the solid materials plays an important role in the classification
of maters and the understanding of physical properties. Then, many topological
materials, such as, topological insulators Kane and Mele (2005); Bernevig and
Zhang (2006); Fu (2009), Dirac semimetals Liu  _et al._ (2014); Wang  _et al._
(2013), Weyl semimetals Shekhar  _et al._ (2015); Lv  _et al._ (2015); Du  _et
al._ (2016) and nodal-line semimetals An  _et al._ (2019); Hu  _et al._
(2016a); Emmanouilidou  _et al._ (2017); Chen  _et al._ (2020) have been
proposed theoretically and confirmed experimentally. Most of them are free of
strong correlation effects. In the presence of strong electron interactions,
very fruitful topological phases can be expected, such as the topological Mott
Pesin and Balents (2010) or Kondo insulators Dzero  _et al._ (2010, 2012); Lu
_et al._ (2013), topological superconductors Qi and Zhang (2011), and
fractional topological insulators Sheng  _et al._ (2011); Regnault and
Bernevig (2011); Neupert  _et al._ (2011). In order to pursue these exotic
phases, searching for suitable compounds, which are strongly correlated
(usually in $d$ and $f$ orbital systems) and topologically nontrivial,
attracts much attention. For instance, SmB6, as a typical rare-earth mixed
valence compound, in which, has been proposed theoretically as a topological
Kondo insulator Dzero  _et al._ (2010, 2012); Lu  _et al._ (2013) and recently
has been confirmed by transport Wolgast  _et al._ (2013); Kim  _et al._
(2013); Li  _et al._ (2014); Thomas  _et al._ (2016), photoemission Xu  _et
al._ (2013); Jiang  _et al._ (2013); Denlinger  _et al._ (2014) and scanning
tunneling microscope (STM) Yee  _et al._ (2013) experiments. At sufficiently
low temperature, the hybridization between $4f$ orbitals and highly dispersive
$5d$ bands results in the formation of “heavy fermion” bands, thus, SmB6 is a
correlated $\mathcal{Z}_{2}$ topological insulator. Recently, a number of
groups Q. S. Wu and Yazyev (2019); Vergniory  _et al._ (2019); Zhang  _et al._
(2019); Tang  _et al._ (2019) predicted thousands of candidate topological
materials by performing systematic high-throughput computational screening
across the databases of known materials. Following these predictions, it is
possible to search for the topological materials in the known compounds
containing $3d$ elements based on the calculations of band structure, and the
measurements of physical properties.

VAs2 crystalizes in a monoclinic structure with space group C12/m1 (No. 12),
as shown in Fig. 1(a), has an isostructure to the transition metal dipnictides
XPn2 (X = Nb, Ta; Pn = P, As, Sb) Xu  _et al._ (2016); Shen  _et al._ (2016);
Wu  _et al._ (2016); Wang  _et al._ (2016a); Li  _et al._ (2016); Yuan  _et
al._ (2016); Luo  _et al._ (2016), which have been widely studied
theoretically and experimentally as a class of topological materials. For
example, Baokai Wang $et$ $al.$ Wang  _et al._ (2019) identified the presence
of a rotational-symmetry-protected topological crystalline insulator (TCI)
states in these compounds based on first-principles calculations combined with
a symmetry analysis. It was found that all these compounds exhibit high
mobilities and extremely large positive MR Shen  _et al._ (2016); Wu  _et al._
(2016); Wang  _et al._ (2016a); Li  _et al._ (2016); Yuan  _et al._ (2016);
Luo  _et al._ (2016). Interestingly, a negative longitudinal MR when the
applied field is parallel to the current direction was observed in both TaSb2
Li  _et al._ (2016) and TaAs2 Luo  _et al._ (2016), similar to that observed
in the known weyl semimetals TaAs family Lv  _et al._ (2015); Wang  _et al._
(2016b); Shekhar  _et al._ (2015); Sun  _et al._ (2015); Lee  _et al._ (2015);
Ghimire  _et al._ (2015); Du  _et al._ (2016); Hu  _et al._ (2016b); Xu  _et
al._ (2015). Compared to the $4d$/$5d$ electrons in NbAs2/TaAs2, the more
localization of the $3d$ electrons in VAs2 might lead to a stronger electronic
correlation, or introduce an additional magnetic scattering to the carriers
(Kondo effect).

In this paper, we grew successfully VAs2 crystals with a monoclinic structure
and measured its longitudinal resistivity $\rho_{xx}(T,H)$, Hall resistivity
$\rho_{xy}(T,H)$, and magnetic susceptibility as a function of temperature at
various magnetic fields, as well as calculated its electronic band structure
and Fermi surface (FS). The band structure calculations show that VAs2 is also
a nodal-line semimetal when spin-orbit coupling (SOC) is ignored. It was found
that the $\rho_{xx}(T)$ measured at $H$ = 0 exhibits a minimum at around 11 K,
which is considered of originate from the magnetic impurities V4+ ($S$ = 1/2)
scattering ($i.e.$ Kondo effect) , evidenced by both the susceptibility
measurements and the fitting of $\rho_{xx}(T)$ data at low temperatures. We
further reveal a nearly quadratic field dependence of MR, reaching 649% at 10
K and 9 T, a field-induced up-turn behavior of $\rho_{xx}(T)$ in this
material, which are attributed to a perfect charge-carrier compensation,
evidenced by both the calculations relying on the FS topology and the Hall
resistivity measurements.

Figure 1: (a) Crystal structure of VAs2. (b) Photograph of a VAs2 crystal with 4 wires for resistance measurements. (c) Single crystal XRD pattern of VAs2. (d) The calculated Fermi surface of VAs2. (e) The Brillouin zone. (f) and (g) Band structures calculated without and with considering SOC. (h) and (i) Nodal lines in the first Brillouin zone from different perspective.  Figure 2:  (a) Temperature dependence of resistivity $\rho(T)$, inset: the $\rho(T)$ data below 40 K and the fitting by Eq. (2). (b) magnetic susceptibility, measured at 1 T, inset: 1/$\chi$(_T_).  Figure 3: (Color online) (a) Temperature dependence of resistivity measured at various magnetic fields. (b) The MR vs. temperature under various magnetic fields. The inset is normalized MR. (c) Temperature dependence of resistivity at 0 and 9 T and their difference. The red and blue lines are the fitting lines using the Kohler scaling law.  Figure 4:  (a) Field dependence of MR at various temperature. (b) MR plotted as a log scale as a function of _H_ /$\rho$xx(0).  Table 1: The obtained parameters by the fitting to $\rho(T)$ data using Eq. (2). $\rho_{0}$ ($\mu\Omega$ cm) |  $\rho_{b}$ ($\mu\Omega$ cm) |  $b$ ($\mu\Omega$ cm/K2) |  $c$ ($\mu\Omega$ cm/K5) |  $T_{K}$ (K) |  $T_{W}$ (K) | $S$
---|---|---|---|---|---|---
2.33 | 4.7 | 1.3$\times$10-3 | 1.0$\times$10-9 | 7.23 | 0.78 | 0.5
Figure 5:  (a) Field dependence of Hall resistivity $\rho_{xy}$ measured at
various temperatures for VAs2 crystal. (b) Charge-carrier mobilities $\mu_{e}$
and $\mu_{h}$, and (inset) the carrier concentrations, $n_{e}$ and $n_{h}$, as
a function of temperature extracted from the two-carrier model analysis of
both $\sigma_{xy}$ and $\sigma_{xx}$ data. Components of the conductivity
tenser, $i.e.$ $\sigma_{xy}$ and $\sigma_{xx}$ in the panels (c) and (d),
respectively, as a function of magnetic field for different temperatures ($<$
60 K). Hollow dots represent experimental data and solid lines are the fitting
curves by using the two-carrier model.

VAs2 single crystals were grown by a chemical vapor transport method. High
purity V and As powder were mixed in a mole ratio 1 : 2, then sealed in an
evacuated silica tube containing I2 as a transport agent with 10 mg/cm3. The
quartz tube was placed in a tube furnace with a temperature gradient of 950∘C
- 750∘C and heated for two weeks. Polyhedral crystals were obtained at the
cold end of the tube. A single crystal with a dimension of $2.0\times
0.73\times 0.29$ mm3 [see Fig. 1(b)] and a cleavage surface (001) was selected
for the transport and magnetic property measurements. The composition was
detected to be V : As = 32.7 : 67, using the Energy Dispersive X-ray
Spectrometer (EDXS), the crystal structure was determined by the single-
crystal X-Ray diffraction, as shown in Fig. 1(c), from which, the lattice
parameter, $c$ = 7.481(7) Å, consistent with the result in Ref. Hulliger
(1964). The longitudinal resistivity and Hall resistivity measurements were
carried out on a physical property measurement system (Quantum Design, PPMS-9
T) with a standard four-probe method [see Fig. 1(b)]. The magnetization
measurements were carried out on a magnetic property measurement system
(Quantum Design, MPMS-7 T). The band structure was calculated by using the
Vienna $ab$ $initio$ simulation package (VASP) Kresse and Furthmüller (1996);
Kresse and Joubert (1999) with a generalized gradient approximation (GGA) of
Perdew, Burke and Ernzerhof (PBE) Perdew  _et al._ (1997) for the exchange
correlation potential. A cutoff energy of 520 eV and a $10\times 10\times 6$
$k$-point mesh were used to perform the bulk calculations. The nodal-line
search and FS calculations were performed by using the open-source software
WannierTools Wu  _et al._ (2018) that is based on the Wannier tight-binding
model (WTBM) constructed using Wannier90 Mostofi  _et al._ (2014).

First at all, we discuss the results of our electronic structure calculations
that extend the initial prediction of the high symmetry line semimetal for
VAs2 Zhang  _et al._ (2019). In order to address the topological character of
VAs2, we calculate its band structure and the FS. As shown in Fig. 1(d), there
are six different FS sheets: two hole-like surfaces (green) at the L and Y
points, and four electron-like (red) surfaces near L and Z points, the volume
of the electron and hole pockets is roughly the same, indicating that VAs2 is
nearly an electron-hole compensated semimetal, also evidenced by the Hall
resistivity measurements discussed as follows. The bulk band structure of VAs2
without and with SOC is presented in Fig. 1(f) and 1(g), respectively. It can
be seen clearly that the bands near the Fermi level mainly arise from the $d$
orbits of V atoms, and the valence and conduction bands cross along the Y-X1,
Z-I1 and L-I high symmetry directions. Without SOC, in the Brillouin zones
(BZ), the nodal lines can be found by using the open-source software
WannierTool Wu  _et al._ (2018), see in Fig. 1(h). There are two type of nodal
lines, one is two nonclosed spiral lines extending across the BZ through point
Z, another is two nodal loops near the L point. When SOC is included, these
nodal lines are gapped [see Fig. 1(g)] and lead to a band inversion along the
Y-X1, Z-I1 and L-I high symmetry directions, driving into a topological
crystalline insulator (TCI), as discussed by Baokai Wang $et$ $al.$ Wang  _et
al._ (2019) for the identical structure transition metal dipnitides RX2 (R =
Nb or Ta; X = P, As or Sb). Although the opening of a local band gap between
the valence and conduction bands occurs when SOC is included, VAs2 preserves
its semimetal character with the presence of electron and hole pockets,
similar to that in NbAs2 reported in Ref. Wang  _et al._ (2019).

Next, we focus on the Kondo effect emerging in the longitudinal resistivity
for VAs2. Figure 2(a) shows the temperature dependence of resistivity,
$\rho_{xx}(T)$. With decreasing temperature, the resistivity, $\rho_{xx}$,
decreases monotonously from $\rho_{xx}$(300 K) = 78 $\mu\Omega$ cm, reaches a
minimum at around 11 K, then increases a little to 6.4 $\mu\Omega$ cm at 2 K,
thus the residual resistivity ratio (RRR) $\rho_{xx}$(300 K)/$\rho_{xx}$ (2 K)
$\sim 12$, much smaller than that observed in the other transition metal
dipnictide crystals, such as NbAs2 ($\sim 75$) Wang  _et al._ (2016a), TaAs2
($\sim 100$) Luo  _et al._ (2016). As we know, the system, containing magnetic
impurities which usually scatter the conduction electrons through the $s-d$
exchange interactions, exhibits a minimum in resistivity at lower
temperatures, $i.e.$ termed as the Kondo effect Kästner and Wassermann (1977);
Liang and Tsuei (1973). As discussed by S. Barua $et$ $al.$ for VSe2 Barua
_et al._ (2017), considering the correction to the resistivity due to the
Ruderman-Kittel-Kasuya-Yosida (RKKY) interactions between the paramagnetic V4+
ions, the Kondo resistivity described by the Hamann expression is modified to:

$\displaystyle\rho_{sd}=\frac{\rho_{0}}{2}\left[1-\frac{\ln(T_{eff}/T_{K})}{\left[\ln^{2}(T_{eff}/T_{K})+S(S+1)\pi^{2}\right]^{1/2}}\right]$
(1)

where $\rho_{0}$ is the unitarity limit, $T_{K}$ is the Kondo temperature and
$S$ is the spin of magnetic impurity, the effective temperature
$T_{eff}=(T^{2}+T^{2}_{W})^{1/2}$, in which $k_{B}T_{W}$ is the effective RKKY
interaction strength Kästner and Wassermann (1977), replaces to $T$ in the
original expression Barua  _et al._ (2017). The temperature dependence of
resistivity, $\rho_{xx}(T)$, is expressed as :

$\displaystyle\rho(T)=\rho_{sd}+bT^{2}+cT^{5}+\rho_{b}$ (2)

where _bT_ n term is the electron-phonon scattering contribution, $\rho_{b}$
is an independent resistivity. We used Eq. (2) to fit the resistivity data
measured at low temperatures (2 - 40 K). The results are shown in the inset in
Fig. 2(a), it is clear that Eq. (2) can well describe the $\rho(T)$ data below
40 K, and the obtained parameters in Eq. (1) and Eq. (2) from the fitting are
listed in TABLE I. The existence of V4+ ($S$ = 1/2) impurities in our VAs2
crystal was confirmed by the magnetic susceptibility measurements. The
temperature dependence of magnetic susceptibility, $\chi(T)$, measured at 1 T
with a field-cooling (FC) process is presented in Fig. 2(b). With decreasing
temperature, the susceptibility $\chi$ decreases a little, reaches a minimum
at around 180 K, then increases strikingly below 100 K. No magnetic transition
was observed in the whole temperature range (2 - 300 K), and the significant
increase of $\chi$ in the lower temperatures was considered to originate from
the contribution of V4+ ($S$ = 1/2) impurities existing in crystal as the
interstitial ions, as observed in VSe2 crystals Barua  _et al._ (2017). We
used the Curie-Weiss law $\chi=\frac{C}{T-\theta}$, to fit the data below 40
K, as shown in Fig. 2(a), the Curie constant $C$ = 1.27
($\pm$0.02)$\times$10-3 emu K/mol, and the curie temperature $\theta$ = -2.61
($\pm$0.08) K, were obtained. The nearly linear relationship between 1/$\chi$
and $T$ below 40 K [see the inset in Fig. 2(b)] demonstrates the reliability
of the fitting. The V4+ ($S$ = 1/2) impurity molar fraction was estimated to
be of 0.34 ($\pm$0.02) $\%$, small amount impurities existing in the crystals.

Third, we discuss the magnetoresisitance (MR) occurring in the nodal-line
semimetal VAs2, with the presence of Kondo scattering of V4+. Figure 3(a)
presents the temperature dependence of longitudinal resistivity,
$\rho_{xx}(T)$, measured at various magnetic fields $H$, with current $I$
applied in the (001) plane, and $H$ $\perp$ (001) plane. Similar to many other
nontrivial and trivial topological semimetals Shen  _et al._ (2016); Wu  _et
al._ (2016); Li  _et al._ (2016), VAs2 also exhibits a large MR. As shown in
Fig. 3(a), an up-turn in $\rho_{xx}$(T) curves under applied magnetic field
occurs at low temperatures: $\rho_{xx}$ increases with decreasing $T$ and then
saturates. Figure 3(b) shows MR as a function of temperature measured at
various magnetic fields, with the conventional definition
$MR$=$\frac{\Delta\rho}{\rho(0)}$=$[\frac{\rho(H)-\rho(0)}{\rho(0)}]$ $\times$
100$\%$. The normalized MR, shown in the inset of Fig. 3(b), has the same
temperature dependence for various fields, excluding the suggestion of a
field-induced metal-insulator transition Zhao  _et al._ (2015); Khveshchenko
(2001) at low temperatures, as discussed in our works addressing the
topological trivial semimetal $\alpha$-WP2 Du  _et al._ (2018), for the nodal-
line semimetal MoO2 Chen  _et al._ (2020), and the work of Thoutam $et$ $al.$
on the type-II weyl semimetal WTe2 Thoutam  _et al._ (2015). Figure 3(c)
displays $\rho_{xx}(T)$, measured at 0 and 9 T, as well as their difference
$\Delta$$\rho$xx = $\rho$xx(_T_ , 9 T) - $\rho$xx(_T_ , 0 T). It is clear that
the resistivity in a applied magnetic field consists of two components,
$\rho_{0}(T)$ and $\Delta\rho_{xx}$, with opposite temperature dependencies.
As discussed by us for $\alpha$-WP2 Du  _et al._ (2018), MoO2 Chen  _et al._
(2020) and by Thoutam $et$ $al$. for WTe2 Thoutam  _et al._ (2015), the
resistivity can be written as :

$\displaystyle\rho_{xx}(T,H)=\rho_{0}(T)[1+\alpha(H/\rho_{0})^{m}]$ (3)

The second term is the magnetic-field-induced resistivity $\Delta\rho_{xx}$,
which follows the Kohler rule with two constants $\alpha$ and $m$.
$\Delta\rho_{xx}$ is proportional to 1/$\rho_{0}$ (when $m$ = 2) and competes
with the first term upon changing temperature, possibly giving rise to a
minimum in $\rho(T)$ curves. Figure 4(a) presents MR as a function of field at
various temperatures. The measured MR is large at low temperatures, reaching
649 $\%$ at 10 K and 9 T, and does not show any of saturation up to the
highest field (9 T) in our measurements. As discussed above, MR can be
described by the Kohler scaling law Wang  _et al._ (2015); Ziman (2001):

$\displaystyle\textit{MR}=\frac{\Delta\rho_{xx}(T,H)}{\rho_{0}(T)}=\alpha(\textit{H}/\rho_{0})^{m}$
(4)

As shown in Fig. 4(b), all MR data from 2 - 100 K collapse onto a single
straight line in the plotted as $MR\sim H/\rho_{0}$ curve, and $\alpha$ =
0.038 $(\mu\Omega$ cm/T$)^{1.76}$ and $m$ = 1.76 were obtained by fitting. The
nearly quadratic field dependence of MR observed in this nodal-line semimetal
VAs2 is attributed to the electron-hole compensation, evidenced by FS
calculations mentioned above, as well as the Hall resistivity measurements
discussed below, which is a common characteristics for the most topologically
nontrivial and trivial semimetals Shen  _et al._ (2016); Wu  _et al._ (2016);
Li  _et al._ (2016). Figure 5(a) displays the Hall resistivity,
$\rho_{xy}(H)$, measured at various temperatures for a VAs2 crystal with $H$
$\parallel$ $c$ axis. The nonlinear field dependence of $\rho_{xy}(H)$ below
100 K demonstrates its semimetal characteristics, in which both electron and
hole carriers coexist. Following the analysis of Ref. Chen  _et al._ (2016)
for $\gamma-$MoTe2 Zhou  _et al._ (2016) and by us for MoO2 Chen  _et al._
(2020), we analyze the longitudinal and Hall resistivity by using the two-
carrier model. In this model, the conductivity tensor, in its complex
representation, is giving by Ali  _et al._ (2014):

$\displaystyle\sigma=\frac{en_{e}\mu_{e}}{1+i\mu_{e}\mu_{0}H}+\frac{en_{h}\mu_{h}}{1+i\mu_{h}\mu_{0}H}$
(5)

where $n_{e}$ ($n_{h}$) and $\mu_{e}$ ($\mu_{h}$) denote the carrier
concentrations and mobilities of electrons (holes), respectively. To
appropriately evaluate the carrier densities and mobilities, we calculated the
Hall conductivity $\sigma_{xy}$ =
-$\frac{\rho_{xy}}{\rho_{xx}^{2}+\rho_{xy}^{2}}$ and the longitudinal
conductivity $\sigma_{xx}$ = $\frac{\rho_{xx}}{\rho_{xx}^{2}+\rho_{xy}^{2}}$
by using the original experimental $\rho_{xy}(H)$ and $\rho_{xx}(H)$ data.
Then, we fit both $\sigma_{xy}(H)$ and $\sigma_{xx}(H)$ data by using the same
fitting parameters and the field dependence given by Zhou  _et al._ (2016):

$\displaystyle\sigma_{xy}=\frac{e(\mu_{0}H)n_{h}\mu_{h}^{2}}{1+\mu_{h}^{2}(\mu_{0}H)^{2}}-\frac{e(\mu_{0}H)n_{e}\mu_{e}^{2}}{1+\mu_{e}^{2}(\mu_{0}H)^{2}}$
(6)

$\displaystyle\sigma_{xx}=\frac{en_{h}\mu_{h}}{1+\mu_{h}^{2}(\mu_{0}H)^{2}}+\frac{en_{e}\mu_{e}}{1+\mu_{e}^{2}(\mu_{0}H)^{2}}$
(7)

Figures 5(c) and 5(d) display the fitting of both the $\sigma_{xy}(H)$ and
$\sigma_{xx}(H)$ measured at $T$ = 2 - 60 K, respectively. The excellent
agreement between our experimental data and the two-carrier model over a broad
range of temperature, confirms the coexistence of electrons and holes in VAs2.
Figure 5(b) shows the obtained $n_{e}$, $n_{h}$, $\mu_{e}$ and $\mu_{h}$
values by fittings as a function of temperature. The almost same values of
$n_{e}$ and $n_{h}$ below 60 K [see the inset Fig. 5(b)], such as $n_{e}$ =
1.77 $\times$1020 cm-3 and $n_{h}$ = 1.69 $\times$1020 cm-3 at 2 K, indicate
that VAs2 is indeed a electron-hole compensated semimetal, consistent with the
above results from the calculation FS. Both electron and hole densities are
estimated to be 1020 cm-3 in our VAs2 crystal, significantly higher than those
in Dirac semimetals, such as, Cd3As2 ($\sim$ 1018 cm-3 Neupane  _et al._
(2014)) and Na3Bi ($\sim$ 1017 cm-3 Xiong  _et al._ (2015)), but comparable to
these of another nodal-line semimetals ZrSiS ($\sim$ 1020 cm-3 Hu  _et al._
(2016a)), and MoO2 ($\sim$ 1020 cm-3 Chen  _et al._ (2020)), demonstrating
further the nodal-line characteristics of VAs2. As shown in Fig. 5(b), it is
clear that the hole mobility $\mu_{h}$ is higher than $\mu_{e}$ in the whole
temperature range (2 - 300 K), such as, at 2 K, $\mu_{h}$ = 4.08 $\times$103
cm2/Vs, $\mu_{e}$ = 1.54 $\times$103 cm2/Vs, but smaller one order of
magnitude than these observed by us in the nodal-line semimetal MoO2 ($\sim$
104 cm2/Vs) Chen  _et al._ (2020), and both $\mu_{h}$ and $\mu_{e}$ decrease
notably with increasing temperature due to the existence of phonon thermal
scattering at higher temperatures. It is worth noting that both $\mu_{h}$ and
$\mu_{e}$ have a little decrease below 11 K, corresponding to the Kondo
scattering from V4+ magnetic impurities mentioned above.

In summary, we calculated the electronic structure, and measured the
longitudinal resistivity, Hall resistivity and magnetic susceptibility for
VAs2. It was found that VAs2 exhibits many common characteristics of the
topological nontrivial/trivial semimetals, such as a large MR reaching 649$\%$
at 10 K and 9 T, a nearly quadratic field dependence of MR, and a field-
induced up-turn behaviour in $\rho$xx(_T_). Both the FS calculations and the
Hall resistivity measurements verify these properties being attributed to a
perfect carrier compensation. Interestingly, the Kondo scattering due to the
existence of V4+ ($S$ = 1/2) magnetic impurities in our VAs2 crystals occurs,
indicating that VAs2 crystal can be used to study the Kondo effect in the
nodal-line semimetal.

ACKNOWLEDGEMENTS: This research is supported by the Ministry of Science and
Technology of China under Grant No. 2016YFA0300402 and the National Natural
Science Foundation of China (NSFC) (NSFC-12074335 and No. 11974095), the
Zhejiang Natural Science Foundation (No. LY16A040012) and the Fundamental
Research Funds for the Central Universities.

## References

  * Laughlin (1983) R. B. Laughlin, Phys. Rev. Lett. 50, 1395 (1983).
  * Kane and Mele (2005) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005).
  * Bernevig and Zhang (2006) B. A. Bernevig and S.-C. Zhang, Phys. Rev. Lett. 96, 106802 (2006).
  * Fu (2009) L. Fu, Phys. Rev. Lett. 103, 266801 (2009).
  * Liu  _et al._ (2014) Z. Liu, B. Zhou, Y. Zhang, Z. Wang, H. Weng, D. Prabhakaran, S.-K. Mo, Z. Shen, Z. Fang, X. Dai, Z. Hussain, and Y. Chen, Science (New York, N.Y.) 343, 864 (2014).
  * Wang  _et al._ (2013) Z. Wang, H. Weng, Q. Wu, X. Dai, and Z. Fang, Phys. Rev. B 88, 125427 (2013).
  * Shekhar  _et al._ (2015) C. Shekhar, A. K. Nayak, Y. Sun, M. Schmidt, M. Nicklas, I. Leermakers, U. Zeitler, Y. Skourski, J. Wosnitza, Z. Liu, _et al._ , Nature Physics 11, 645 (2015).
  * Lv  _et al._ (2015) B. Q. Lv, H. M. Weng, B. B. Fu, X. P. Wang, H. Miao, J. Ma, P. Richard, X. C. Huang, L. X. Zhao, G. F. Chen, Z. Fang, X. Dai, T. Qian, and H. Ding, Phys. Rev. X 5, 031013 (2015).
  * Du  _et al._ (2016) J. Du, H. Wang, Q. Chen, Q. Mao, R. Khan, B. Xu, Y. Zhou, Y. Zhang, J. Yang, B. Chen, _et al._ , Science China Physics, Mechanics & Astronomy 59, 657406 (2016).
  * An  _et al._ (2019) L. An, X. Zhu, W. Gao, M. Wu, W. Ning, and M. Tian, Phys. Rev. B 99, 045143 (2019).
  * Hu  _et al._ (2016a) J. Hu, Z. Tang, J. Liu, X. Liu, Y. Zhu, D. Graf, K. Myhro, S. Tran, C. N. Lau, J. Wei, and Z. Mao, Phys. Rev. Lett. 117, 016602 (2016a).
  * Emmanouilidou  _et al._ (2017) E. Emmanouilidou, B. Shen, X. Deng, T.-R. Chang, A. Shi, G. Kotliar, S.-Y. Xu, and N. Ni, Phys. Rev. B 95, 245113 (2017).
  * Chen  _et al._ (2020) Q. Chen, Z. Lou, S. Zhang, B. Xu, Y. Zhou, H. Chen, S. Chen, J. Du, H. Wang, J. Yang, Q. Wu, O. V. Yazyev, and M. Fang, Phys. Rev. B 102, 165133 (2020).
  * Pesin and Balents (2010) D. Pesin and L. Balents, Nature Physics 6, 376 (2010).
  * Dzero  _et al._ (2010) M. Dzero, K. Sun, V. Galitski, and P. Coleman, Phys. Rev. Lett. 104, 106408 (2010).
  * Dzero  _et al._ (2012) M. Dzero, K. Sun, P. Coleman, and V. Galitski, Phys. Rev. B 85, 045130 (2012).
  * Lu  _et al._ (2013) F. Lu, J. Zhao, H. Weng, Z. Fang, and X. Dai, Phys. Rev. Lett. 110, 096401 (2013).
  * Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Sheng  _et al._ (2011) D. N. Sheng, Z. C. Gu, K. Sun, and L. Sheng, Nature Communications 2, 389 (2011).
  * Regnault and Bernevig (2011) N. Regnault and B. A. Bernevig, Phys. Rev. X 1, 021014 (2011).
  * Neupert  _et al._ (2011) T. Neupert, L. Santos, C. Chamon, and C. Mudry, Phys. Rev. Lett. 106, 236804 (2011).
  * Wolgast  _et al._ (2013) S. Wolgast, i. m. c. b. u. i. e. i. f. Kurdak, K. Sun, J. W. Allen, D.-J. Kim, and Z. Fisk, Phys. Rev. B 88, 180405 (2013).
  * Kim  _et al._ (2013) D. Kim, S. Thomas, T. Grant, J. Botimer, Z. Fisk, and J. Xia, Scientific reports 3, 3150 (2013).
  * Li  _et al._ (2014) G. Li, Z. Xiang, F. Yu, T. Asaba, B. Lawson, P. Cai, C. Tinsman, A. Berkley, S. Wolgast, and Y. S. a. Eo, Science 346, 1208 (2014).
  * Thomas  _et al._ (2016) S. Thomas, D. J. Kim, S. B. Chung, T. Grant, Z. Fisk, and J. Xia, Phys. Rev. B 94, 205114 (2016).
  * Xu  _et al._ (2013) N. Xu, X. Shi, P. K. Biswas, C. E. Matt, R. S. Dhaka, Y. Huang, N. C. Plumb, M. Radović, J. H. Dil, E. Pomjakushina, K. Conder, A. Amato, Z. Salman, D. M. Paul, J. Mesot, H. Ding, and M. Shi, Phys. Rev. B 88, 121102 (2013).
  * Jiang  _et al._ (2013) J. Jiang, S. Li, T. Zhang, Z. Sun, F. Chen, Z. R. Ye, M. Xu, Q. Q. Ge, S. Y. Tan, and X. H. a. Niu, Nature Communications 4, 3010 (2013).
  * Denlinger  _et al._ (2014) J. D. Denlinger, J. W. Allen, J.-S. Kang, K. Sun, B.-I. Min, D.-J. Kim, and Z. Fisk, in  _Proceedings of the International Conference on Strongly Correlated Electron Systems (SCES2013)_ (2014) p. 017038.
  * Yee  _et al._ (2013) M. M. Yee, Y. He, A. Soumyanarayanan, D.-J. Kim, Z. Fisk, and J. E. Hoffman, arXiv preprint arXiv:1308.1085 (2013).
  * Q. S. Wu and Yazyev (2019) N. M. Q. S. Wu, G. Autès and O. V. Yazyev, Materials Cloud Archive 2019.0019/v2 (2019), 10.24435/materialscloud:2019.0019/v2.
  * Vergniory  _et al._ (2019) M. Vergniory, L. Elcoro, C. Felser, N. Regnault, B. A. Bernevig, and Z. Wang, Nature 566, 480 (2019).
  * Zhang  _et al._ (2019) T. Zhang, Y. Jiang, Z. Song, H. Huang, Y. He, Z. Fang, H. Weng, and C. Fang, Nature 566, 475 (2019).
  * Tang  _et al._ (2019) F. Tang, H. C. Po, A. Vishwanath, and X. Wan, Nature 566, 486 (2019).
  * Xu  _et al._ (2016) C. Xu, J. Chen, G.-X. Zhi, Y. Li, J. Dai, and C. Cao, Phys. Rev. B 93, 195106 (2016).
  * Shen  _et al._ (2016) B. Shen, X. Deng, G. Kotliar, and N. Ni, Phys. Rev. B 93, 195119 (2016).
  * Wu  _et al._ (2016) D. Wu, J. Liao, W. Yi, X. Wang, P. Li, H. Weng, Y. Shi, Y. Li, J. Luo, X. Dai, and Z. Fang, Appl. Phys. Lett. 108, 042105 (2016).
  * Wang  _et al._ (2016a) Y.-Y. Wang, Q.-H. Yu, P.-J. Guo, K. Liu, and T.-L. Xia, Phys. Rev. B 94, 041103 (2016a).
  * Li  _et al._ (2016) Y. Li, L. Li, J. Wang, T. Wang, X. Xu, C. Xi, C. Cao, and J. Dai, Phys. Rev. B 94, 121115 (2016).
  * Yuan  _et al._ (2016) Z. Yuan, H. Lu, Y. Liu, J. Wang, and S. Jia, Phys. Rev. B 93, 184405 (2016).
  * Luo  _et al._ (2016) Y. Luo, R. D. McDonald, P. F. S. Rosa, B. Scott, N. Wakeham, N. J. Ghimire, E. D. Bauer, J. D. Thompson, and F. Ronning, Sci. Rep. 6, 27294 (2016).
  * Wang  _et al._ (2019) B. Wang, B. Singh, B. Ghosh, W.-C. Chiu, M. M. Hosen, Q. Zhang, L. Ying, M. Neupane, A. Agarwal, H. Lin, and A. Bansil, Phys. Rev. B 100, 205118 (2019).
  * Wang  _et al._ (2016b) Z. Wang, Y. Zheng, Z. Shen, Y. Lu, H. Fang, F. Sheng, Y. Zhou, X. Yang, Y. Li, C. Feng, and Z.-A. Xu, Phys. Rev. B 93, 121112 (2016b).
  * Sun  _et al._ (2015) Y. Sun, S.-C. Wu, and B. Yan, Phys. Rev. B 92, 115428 (2015).
  * Lee  _et al._ (2015) C.-C. Lee, S.-Y. Xu, S.-M. Huang, D. S. Sanchez, I. Belopolski, G. Chang, G. Bian, N. Alidoust, H. Zheng, M. Neupane, B. Wang, A. Bansil, M. Z. Hasan, and H. Lin, Phys. Rev. B 92, 235104 (2015).
  * Ghimire  _et al._ (2015) N. J. Ghimire, Y. Luo, M. Neupane, D. J. Williams, E. D. Bauer, and F. Ronning, Journal of Physics: Condensed Matter 27, 152201 (2015).
  * Hu  _et al._ (2016b) J. Hu, J. Liu, D. Graf, S. Radmanesh, D. Adams, A. Chuang, Y. Wang, I. Chiorescu, J. Wei, L. Spinu, _et al._ , Scientific Reports 6, 18674 (2016b).
  * Xu  _et al._ (2015) S.-Y. Xu, I. Belopolski, N. Alidoust, M. Neupane, G. Bian, C. Zhang, R. Sankar, G. Chang, Z. Yuan, C.-C. Lee, S.-M. Huang, H. Zheng, J. Ma, D. S. Sanchez, B. Wang, A. Bansil, F. Chou, P. P. Shibayev, H. Lin, S. Jia, and M. Z. Hasan, Science 349, 613 (2015).
  * Hulliger (1964) F. Hulliger, Nature 204, 775 (1964).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Perdew  _et al._ (1997) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 78, 1396 (1997).
  * Wu  _et al._ (2018) Q. Wu, S. Zhang, H.-F. Song, M. Troyer, and A. A. Soluyanov, Computer Physics Communications 224, 405 (2018).
  * Mostofi  _et al._ (2014) A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Computer Physics Communications 185, 2309 (2014).
  * Kästner and Wassermann (1977) J. Kästner and E. Wassermann, Journal of Low Temperature Physics 29, 411 (1977).
  * Liang and Tsuei (1973) V. K. C. Liang and C. C. Tsuei, Phys. Rev. B 7, 3215 (1973).
  * Barua  _et al._ (2017) S. Barua, M. C. Hatnean, M. Lees, and G. Balakrishnan, Scientific reports 7, 1 (2017).
  * Zhao  _et al._ (2015) Y. Zhao, H. Liu, J. Yan, W. An, J. Liu, X. Zhang, H. Wang, Y. Liu, H. Jiang, Q. Li, Y. Wang, X.-Z. Li, D. Mandrus, X. C. Xie, M. Pan, and J. Wang, Phys. Rev. B 92, 041104 (2015).
  * Khveshchenko (2001) D. V. Khveshchenko, Phys. Rev. Lett. 87, 206401 (2001).
  * Du  _et al._ (2018) J. Du, Z. Lou, S. Zhang, Y. Zhou, B. Xu, Q. Chen, Y. Tang, S. Chen, H. Chen, Q. Zhu, H. Wang, J. Yang, Q. Wu, O. V. Yazyev, and M. Fang, Phys. Rev. B 97, 245101 (2018).
  * Thoutam  _et al._ (2015) L. R. Thoutam, Y. L. Wang, Z. L. Xiao, S. Das, A. Luican-Mayer, R. Divan, G. W. Crabtree, and W. K. Kwok, Phys. Rev. Lett. 115, 046602 (2015).
  * Wang  _et al._ (2015) Y. L. Wang, L. R. Thoutam, Z. L. Xiao, J. Hu, S. Das, Z. Q. Mao, J. Wei, R. Divan, A. Luican-Mayer, G. W. Crabtree, and W. K. Kwok, Phys. Rev. B 92, 180402 (2015).
  * Ziman (2001) J. Ziman,  _Electrons and Phonons, Classics Series_ (Oxford University Press, 2001).
  * Chen  _et al._ (2016) F. C. Chen, H. Y. Lv, X. Luo, W. J. Lu, Q. L. Pei, G. T. Lin, Y. Y. Han, X. B. Zhu, W. H. Song, and Y. P. Sun, Phys. Rev. B 94, 235154 (2016).
  * Zhou  _et al._ (2016) Q. Zhou, D. Rhodes, Q. R. Zhang, S. Tang, R. Schönemann, and L. Balicas, Phys. Rev. B 94, 121101 (2016).
  * Ali  _et al._ (2014) M. N. Ali, J. Xiong, S. Flynn, J. Tao, Q. D. Gibson, L. M. Schoop, T. Liang, N. Haldolaarachchige, M. Hirschberger, N. P. Ong, _et al._ , Nature 514, 205 (2014).
  * Neupane  _et al._ (2014) M. Neupane, S.-Y. Xu, R. Sankar, N. Alidoust, G. Bian, C. Liu, I. Belopolski, T.-R. Chang, H.-T. Jeng, H. Lin, _et al._ , Nature communications 5, 1 (2014).
  * Xiong  _et al._ (2015) J. Xiong, S. K. Kushwaha, T. Liang, J. W. Krizan, M. Hirschberger, W. Wang, R. J. Cava, and N. P. Ong, Science 350, 413 (2015).