††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

# Stress-Induced Mutagenesis
Can Further Boost Population Success in Static Ecology

Kien T. Pham  Department of Aerospace Engineering, School of Transportation
Engineering, Hanoi University of Science and Technology, 1 Dai Co Viet Str.,
Hai Ba Trung District, Hanoi 100000, Vietnam.  Duc M. Nguyen  The College,
University of Chicago, Chicago, IL 60637, USA.  Duy V. Tran  University of
Technology - VNUHCM, 228 Ly Thuong Kiet, Ho Chi Minh 700000, Vietnam.  Vi D.
Ao  University of Science - VNUHCM, 227-Nguyen Van Cu Street, 5th District, Ho
Chi Minh City, Vietnam.  Huy D. Tran  Department of Physics, The Hong Kong
University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong,
P.R. China.  Tuan K. Do  Department of Mathematics, UCLA, Los Angeles, CA
90095-1555, USA.  Trung V. Phan<EMAIL_ADDRESS>Department of Molecular,
Cellular, and Developmental Biology, Yale University, New Haven, CT 06520, USA

###### Abstract

We have developed a mathematical model that captures stress-induced
mutagenesis, a fundamental aspect of pathogenic and neoplastic evolutionary
dynamics, on the fitness landscape with multiple relevant genetic traits as a
$\mathscrsfs{D}$-dimensional Euclidean space. In this framework, stress-
induced mutagenesis manifests as a heterogeneous diffusion process. We show
how increasing mutations, and thus reducing exploitation, in a static ecology
with fixed carrying capacity and maximum growth rates, can paradoxically boost
population size. Remarkably, this unexpected biophysical phenomenon applies
universally to any number of traits $\mathscrsfs{D}=1,2,3,...$

##  I Introduction

The concepts of exploration and exploitation are often thought of as opposing
strategies, with finding the right balance between the two being the key to
achieving success, both in the natural world and in human endeavors Gupta
(2006); Uotila (2009); Berger-Tal (2014); Guisado-González (2017); Alba
(2005). Exploitation focuses on maximizing the utilization of existing
resources and knowledge to achieve optimal outcomes. In contrast, exploration
is the process of searching for better options, which becomes crucial when
facing ever-changing challenges. On the other hand, if the challenge remains
static, and the optimal solution is well-defined globally, then exploration
may be perceived as futile. In biology, genetic mutation functions as the
exploration mechanism in the extensive range of all possible combinations of
genetic traits (the biological landscape) Wright (1932); Kimura (1962, 1964),
and a species’ success can be gauged by how large its population size is Phan
(2021). For an unchanging environment, it is suggested that for maximum
success, the population distribution in the landscape should be as localized
as possible around the genetic traits with the highest growth rate, rendering
a lower mutation rate more desirable.

However, this claim, as we find out, is fundamentally flawed.

In natural complex and dynamic ecology, where living systems strongly interact
with environments Phan (2020), there exists a coupling between the abstract
space of genetic trait variations and outside physical space. We reveal that
this coupling can have a significant impact through stress-induced
mutagenesis, a common evolutionary mechanism that has been observed not only
at the microbial level Bjedov (2003) but also in cancer cells Fitzgerald
(2017). This curious mechanism is triggered by high biological stress such as
starvation or exposure to antibiotics Imasheva (1999); Hoffmann (2000). For
E.coli bacteria, this process involves the engagement of the SOS response,
leading to the induction of low-fidelity error-prone replication polymerases
and a sharp increased mutation rate during DNA replication from the normally
low value of $D_{l}\propto 10^{-9}$ to a high rate of $D_{h}\propto 10^{-5}$
mutations per base pair per generation Cirz (2005); Bos (2015). In general,
there are many mechanisms by which genetic change can be produced when
organisms are under stress Foster (2007). The importance of stress-induced
mutagenesis, as an essential evolutionary dynamic in the realm of
microbiology, has been demonstrated by laboratory studies showing that at
least $80\%$ of natural isolates of E.coli from diverse environments worldwide
can exhibit this behavior Pribis (2022). For humanity, it poses a significant
threat e.g. the direct contribution to the rapid development of drug
resistance Zhang (2011); Wu (2014); Li (2021), leading to complications in
treating infectious diseases. According to the World Health Organization, an
estimated 1.27 million deaths worldwide resulted from antibiotic-resistant
infections in 2019 Murray (2022).

Previously, we have demonstrated that more exploration non-stress-induced can
lead to better population success for living system in dynamic surrounding
Phan (2021). It also has been reported in the literature that stress-induced
mutagenesis can drastically enhance adaptation to ever-changing environment
Foster (2007); Ram (2014). In this work, we find that the benefit of stress-
induced mutagenesis is even more fundamental: the population success $S=B/K$
Phan (2021) can already be boosted in a static ecology where the carrying
capacity $K$ Tsoularis (2002) and the maximum growth rates $R(\vec{x})$
Gavrilets (2018); Wright (1932) for all combinations $\vec{x}$ of
$\mathscrsfs{D}$ relevant genetic traits (e.g. scalar-strategies Vincent
(2005)) on the Euclidean-topology landscape
($\vec{x}\in\mathbb{R}^{\mathscrsfs{D}}$) are constants with time. $B=\int
d^{\mathscrsfs{D}}\vec{x}b(\vec{x})$, the integration of population density
distribution $b(\vec{x})$ over all this abstract space, is the total
population in the ecology. Here we present a simple mathematical model in
which stress-induced mutagenesis is described on the landscape by a diffusion
process Wright (1932); Kimura (1962, 1964) with growth-dependent diffusivity
$D[G]$, where the growth is the same as in the simple logistic-growth model
$G=R(1-S)$ Tsoularis (2002) and upping mutation rate happens sharply at no-
growth $G=0$ (e.g. for E.coli, this is when cell stops dividing and gets
elongated Bos (2015); Phan (2018)):

$D[G>0]=D_{l}\ ,\ D[G<0]=D_{h}\geq D_{l}\ \ ,$ (1)

analogous to the realistic response Cirz (2005); Bos (2015). Let us define:

$\epsilon=(D_{h}-D_{l})/D_{l}\geq 0\ ,$ (2)

then $\epsilon$ can be used as a metric to quantify the strength of stress-
induced mutagenesis, which increases when this effect becomes stronger. We
show, with analytical and numerical investigations, that the success $S$ at
stationary-state goes up monotonically with the strength $\epsilon$.
Remarkably, we also can prove the universality of this finding for any number
of traits $\mathscrsfs{D}\in\mathbb{N}$, indicating the generality of our
claim. Our results highlight an underdeveloped sector of theoretical
evolutionary dynamics, which is common in the world of cells and can exhibit
such strange and counter-intuitive phenomena, and emphasize the potential of
ultilizing this mathematical framework to uncover the emergent complexities of
evolution.

##  II Theoretic Considerations for Contributions from Stress-Induced
Mutagenesis

In Section II.1, we introduce our mathematical toy-model for evolutionary
dynamics, which incorporates genetic diversity with varying fitness, growth
constraints that reflect ecological limitations on resources and physical
competition, and, most importantly, stress-induced mutagenesis. Moving to
Section II.2, we demonstrate that without stress-induced mutagenesis
$\epsilon=0$, increasing the exploration rate by uniformly raising the
mutation rate $D_{l}=D_{h}\uparrow$ leads to a decrease in population success
$S\downarrow$. Finally, in Section II.3, we examine the effect of stress-
induced mutagenesis $\epsilon>0$, which promotes exploration only for unfit
combinations of genetic traits with a high mutation rate $D_{h}\uparrow$. We
reveal how this evolutionary mechanism can substantially and universally boost
population success $S\uparrow$, even rescuing the population from total
extinction.

###  II.1 An Integro-differential equation for Stress-Induced Mutagenesis

Mutation can be regarded as a stochastic process involving random-walk on the
abstract space of genetic variations Kimura (1964); Phan (2021). In this
space, each combination of traits is represented by a point $\vec{x}$
($\vec{x}\in\mathbb{R}{D}$ for $\mathscrsfs{D}$ evolutionary scalar strategies
Vincent (2005)) and exhibits unchanging maximum growth rate $R(\vec{x})$ in
static ecology. The distribution density $b(\vec{x},t)$ of the population in
this landscape can be modeled by the following Fokker-Planck equation Risken
(1996):

$\partial_{t}b=\nabla^{2}\left(Db\right)+Rb\ ,$ (3)

where the diffusivity $D$ represents the speed of mutations.

In every ecological system, unlimited population growth is unsustainable due
to limited resources in the environment. The logistic model of population
growth is a fundamental description that accounts for this limitation
Tsoularis (2002). In such models, competition for resources results in a
carrying capacity $K$ Getz (1991, 1994), which can be incorporated in the
mathematical framework by modifying Eq. (3) into an integro-differential
equation:

$\partial_{t}b=\nabla^{2}\left(Db\right)+Gb\ ,$ (4)

where $G=R\left(1-S\right)$ is the growth rate, $S=B/K$ is the population
success, and $B$ is the population size which can be calculated as the
integration of population density distribution $B(t)=\int
d^{\mathscrsfs{D}}\vec{x}b(\vec{x},t)$. To include a realistic representation
of stress-induced mutagenesis, we use the two-state growth-dependent
diffusivity $D[G]$ as argued in (1). To increase the veracity of the model, it
is necessary to apply a generalized Richards law Richards (1959) to regulate
growth, where the effective exponent varies with population size Swartz
(2022); Tran (2022).

Here comes the “toy” parts of our toy-model. For the sake of simplicity we
assume that the fitness landscape is symmetric around the position of the
optimal combinations of traits (taken to be the origin) $\vec{x}_{op}=0$:

$R(\vec{x})\approx
R_{0}\left[1-\left(\frac{\vec{x}}{\lambda}\right)^{2}\right]\ .$ (5)

The further away from the optimal combinations, the lower the fitness.
Negative growth means the number of deaths exceeds the number of births, a
crucial indicator on the fitness landscape showing the impact of natural
selection McCandlish (2011). Since $G\propto R$, the characteristic sharp
transition from $D_{l}$ to $D_{h}$ in $D[G]$ at no-growth also happens right
at $R=0$, where $|\vec{x}|=\lambda$. We call $|\vec{x}|<\lambda$ the fit
region where $R>0$, and $|\vec{x}|>\lambda$ the unfit region where $R<0$ Phan
(2021). The basics of our model are conveyed in Fig. 1.

Figure 1: The landscape of our toy-model, illustrated for a single genetic
trait $\mathscrsfs{D}=1$. The effective diffusivity $D[G]$ shifts from a low
value $D_{l}$ to a high value $D_{h}$ when the growth rate $G\propto R$
switches sign.

We can change the spatial coordinates from Cartesian $\vec{x}$ to generalized-
spherical coordinates $\left(r,\Omega_{\mathscrsfs{D}-1}\right)$ then
subsequently apply a rescaling $\rho=r/\lambda$, so that the dimensionless
radial measure $\rho<1$ stands for the fit region and $\rho>1$ corresponds to
the unfit region. We can also rescale the population measure
$\beta=b\lambda^{\mathscrsfs{D}}/K$, so that
$S=\int\rho^{\mathscrsfs{D}-1}d\rho d\Omega_{\mathscrsfs{D}-1}\beta$. Let us
consider a radial symmetric distribution $\beta=\beta(\rho,t)$, Eq. (4)
becomes:

$\partial_{t}\beta=\frac{1}{\rho^{\mathscrsfs{D}-1}}\partial_{\rho}\Bigg{[}\rho^{\mathscrsfs{D}-1}\partial_{\rho}\left(\Big{[}1+\epsilon\Theta(\rho-1)\Big{]}\beta\right)\Bigg{]}+\Xi\left(1-S\right)\left(1-\rho^{2}\right)\beta\
,$ (6)

where $\Theta(\zeta)$ is the Heaviside function in which $\Theta(\zeta<0)=0$
and $\Theta(\zeta>0)=1$ Heaviside (1893), and $\Xi=R_{0}\lambda^{2}/D_{l}$.
The larger fit region size $\lambda\uparrow$, the smaller the parameter
$\Xi\downarrow$. Same goes for higher maximum optimal growth rate
$R_{0}\uparrow$ and lower diffusivity $D_{l}\downarrow$.

Through the solution $\beta=\beta_{st}(\rho)$ for the stationary state of Eq.
(6) by setting $\partial_{t}\beta=0$ and the subsequent determination of the
resulting $S_{st}$, we can gain new insights into the effect of stress-induced
mutagenesis at population level.

###  II.2 Without Stress-Induced Mutagenesis

When there is no stress-induced mutagenesis, $\epsilon=0$ and Eq. (6)
simplifies:

$0=\frac{1}{\rho^{\mathscrsfs{D}-1}}\partial_{\rho}\left(\rho^{\mathscrsfs{D}-1}\partial_{\rho}\beta_{st}\right)+\omega_{l}^{2}\left(1-\rho^{2}\right)\beta_{st},$
(7)

in which we define:

$\omega_{l}^{2}=\Xi(1-S_{st})=\frac{R_{0}\lambda^{2}}{D_{l}}(1-S_{st})\ .$ (8)

The general solution for Eq. (7) are given by the linear superposition of a
Tricomi confluent hypergeometric-function $\operatorname{U}(...)$ and a
Laguerre polynomial $\operatorname{L}(...)$Abramowitz (1948):

$\beta_{st}=e^{-\frac{1}{2}\omega_{l}\rho^{2}}\left[\beta_{1}\
\operatorname{U}\left(\frac{\mathscrsfs{D}-\omega_{l}}{4},\frac{\mathscrsfs{D}}{2},\omega_{l}\rho^{2}\right)\
+\ \beta_{2}\
\operatorname{L}\left(-\frac{\mathscrsfs{D}-\omega_{l}}{4},\frac{\mathscrsfs{D}}{2}-1,\omega_{l}\rho^{2}\right)\right]\
.$ (9)

The conditions that a physical $\beta_{st}$ has to satisfy are: (i) has no
kink at $\rho=0$; (ii) asymptotically goes to $0$ at far-away infinity
$\rho\rightarrow 0$; and (iii) can never be negative anywhere. These
requirements uniquely determine the shape of $\beta_{st}$ and the value of
$\omega_{l}$:

$\beta_{st}\propto e^{-\frac{1}{2}\omega_{l}\rho^{2}}\ ,\
\omega_{l}=\mathscrsfs{D}\ .$ (10)

Therefore, the stationary population success is given by:

$S_{st}=1-\frac{\omega_{l}^{2}}{\Xi}=1-\frac{D_{l}\mathscrsfs{D}^{2}}{R_{0}\lambda^{2}}\
,$ (11)

which decreases as the diffusivity $D_{l}\uparrow$ increases. This finding
corresponds to the expectation that more exploration results in lower success
for static ecology. We note that it is also the same consequence with when
there are more traits $\mathscrsfs{D}\uparrow$.

We show how this intuitive claim fails in the next Section, once effect of
stress-induced mutagenesis $\epsilon>0$ is considered.

###  II.3 With Stress-Induced Mutagenesis

Figure 2: The surface plot $\omega_{l}(\mathscrsfs{D},\epsilon)$, in which the
value of $\omega_{l}$ always decrease as $\epsilon$ increases for every
natural dimensionality $\mathscrsfs{D}$. The parameter space of investigation
spans from $\epsilon=10^{-2}$, $\mathscrsfs{D}=10^{0}$ to $\epsilon=10^{8}$,
$\mathscrsfs{D}=10^{1}$ (see Appendix A).

To deal with Eq. (6) when stress-induced mutagenesis is in effected
$\epsilon>0$, we have to stitch the two solutions $\beta_{st}^{(l)}$ and
$\beta_{st}^{(h)}$ of the form given in Eq. (9) together, in which
$\beta_{st}^{(l)}$ is valid for the fit region $\rho<1$ and $\beta_{st}^{(h)}$
is applicable for the unfit region $\rho>1$. Because in general
$\operatorname{U}(...)$ has a kink at $\rho=0$ and $\operatorname{L}(...)$
diverges when $\rho\rightarrow\infty$, we are looking for a physical
$\beta_{st}$ of the following ansatz:

$\begin{split}\beta_{st}\Big{|}_{\rho<1}&=\beta^{(l)}_{st}\propto
e^{-\frac{1}{2}\omega_{l}\rho^{2}}\
\operatorname{L}\left(-\frac{\mathscrsfs{D}-\omega_{l}}{4},\frac{\mathscrsfs{D}}{2}-1,\omega_{l}\rho^{2}\right)\
,\\\ \beta_{st}\Big{|}_{\rho>1}&=\beta^{(h)}_{st}\propto
e^{-\frac{1}{2}\omega_{h}\rho^{2}}\
\operatorname{U}\left(\frac{\mathscrsfs{D}-\omega_{h}}{4},\frac{\mathscrsfs{D}}{2},\omega_{h}\rho^{2}\right)\
.\end{split}$ (12)

where we define $\omega_{h}$ as a function of $\omega_{l}$, using Eq. (2):

$\omega_{h}^{2}=\left(\frac{\omega_{l}}{\sqrt{1+\epsilon}}\right)^{2}=\frac{R_{0}\lambda^{2}}{D_{h}}(1-S_{st})\
.$ (13)

For regular flux-continuation $\nabla\left(D\beta\right)$ at the transition
interface $\rho=1$, we need to match up Eq. (12) such that:

$D_{l}\beta_{st}\Big{|}_{\rho=1^{-}}=D_{h}\beta_{st}\Big{|}_{\rho=1^{+}}\ ,\
D_{l}\partial_{\rho}\beta_{st}\Big{|}_{\rho=1^{-}}=D_{h}\partial_{\rho}\beta_{st}\Big{|}_{\rho=1^{+}}\
.$ (14)

After performing some algebraic manipulations, it follows that the value of
$\omega_{l}$ must satisfy:

$\omega_{l}\left[1+\frac{2\
\operatorname{L}\left(-\frac{\mathscrsfs{D}-\omega_{l}}{4}-1,\frac{\mathscrsfs{D}}{2},\omega_{l}\right)}{\operatorname{L}\left(-\frac{\mathscrsfs{D}-\omega_{l}}{4},\frac{\mathscrsfs{D}}{2}-1,\omega_{l}\right)}\right]=\omega_{h}\left[1+\frac{\left(\mathscrsfs{D}-\omega_{h}\right)\
\operatorname{U}\left(\frac{\mathscrsfs{D}-\omega_{h}}{4}+1,\frac{\mathscrsfs{D}}{2}+1,\omega_{h}\right)}{2\
\operatorname{U}\left(\frac{\mathscrsfs{D}-\omega_{h}}{4},\frac{\mathscrsfs{D}}{2},\omega_{h}\right)}\right]\
.$ (15)

This equation can be solved numerically, as we describe in Appendix A, giving
us the surface function of $\omega_{l}(\mathscrsfs{D},\epsilon)$ as shown in
Fig. 2:

$\boxed{\omega_{l}(\mathscrsfs{D},0)=\mathscrsfs{D}\ >\
\omega_{l}(\mathscrsfs{D},\epsilon)\Big{|}_{\epsilon>0}\ \
\forall\mathscrsfs{D}\in\mathbb{N}\ .}$ (16)

This is our key result. Since $S_{st}$ and $\omega_{l}$ are monotonically
opposed from the definition in Eq. (8), this means it is universal for any
number of relevant genetic traits $\mathscrsfs{D}=1,2,3,...$ that stress-
induced mutagenesis $\epsilon\uparrow$ always boosts up population success
$S_{st}\uparrow$. For $\mathscrsfs{D}=1$ we get the asymptotic behavior
$\omega_{l}\propto\epsilon^{-1/6}$ at large value of $\epsilon\gg 1$, which is
the regime of stress-induced mutagenesis observed in real living E.coli
bacteria $\epsilon\sim 10^{4}$ Cirz (2005); Bos (2015). In that same regime,
for $\mathscrsfs{D}=2$ the value of $\omega_{l}$ very gradually approaches $0$
and for $\mathscrsfs{D}>2$ the value of $\omega_{l}$ converges to finite
constants. The mathematical details regarding these analytical considerations
can be found in Appendix B.

Although the impact of stress-induced mutagenesis on $S_{st}$ may appear
negligible for biological systems having a wide fit region $\lambda\uparrow$
(associated with a large value of $\Xi\uparrow$, as shown in Fig. 3A for
$\Xi=120$), it is of substantial significance when the fit region is narrow
$\lambda\downarrow$ (resulting in a minuscule value of $\Xi\downarrow$). To
illustrate this point, consider the case of $\Xi=6$, where Fig. 3B displays an
extinction swamp with no possibility for any non-zero stationary population to
thrive Phan (2021); Wang (2022). For $\mathscrsfs{D}=8$, as $\epsilon$
approaching $\sim 10^{4}$ (similar to E.coli Cirz (2005); Bos (2015)) the
population can escape the swamp and become sustainable. In other words,
stress-induced mutagenesis has come to the rescue!

Figure 3: The surface plots for $S(\mathscrsfs{D},\epsilon)$, investigated in
the parameter space spans from $\epsilon=10^{-2}$, $\mathscrsfs{D}=10^{0}$ to
$\epsilon=10^{8}$, $\mathscrsfs{D}=10^{1}$, in which: (A) $\Xi=120$, the
stress-induced mutagenesis effect seems negligible at low number of relevant
genetic traits $\mathscrsfs{D}$. (B) $\Xi=6$, the stress-induced mutagenesis
effect is substantial at low number of traits $\mathscrsfs{D}$. We can observe
an extinction swamp here. We highlight $S(8,\epsilon)$ by a green curve, which
is deep in the swamp until $\epsilon\gtrsim 10^{4}$.

In order to gain insights into the emergence of boosted population success in
response to stress-induced mutagenesis, it is crucial to understand the
distribution of the population across the fitness landscape. To that end, we
show in Fig. 4 the analytic results obtained from equations (12) and (15),
which has also been verified through an agent-based random-walk simulation as
described in detail in Appendix C. Rapid dispersion of the unfit subpopulation
across the fitness landscape alleviates the burden of resource limitation on
fit individuals, enhancing the effective fitness of the population and
increasing the number of fit individuals compared to the absence of stress-
induced mutagenesis.

Figure 4: We compare our analytic results following from Eq. (12) and Eq. (15)
with the stationary states obtained from an agent-based simulation (see
Appendix C). Increasing the strength of stress-induced mutagenesis leads to an
increase in the number of fit individuals and a wider dispersion of the unfit
subpopulation across the landscape.

##  III Discussion

The renowned adage, “nothing makes sense in biology except in the light of
evolution Dobzhansky (1973),” suggests that natural selection favors
evolutionary adaptations that assist in optimizing the success of living
communities within their ecological niche Morris (2017). Nevertheless, the
existence of stress-induced mutagenesis, which may have been widely present at
the dawn of life NRC (1990); Bjedov (2003), poses a baffling question: how
does this up-exploration mechanism enhance population success in seemingly
unchanging environments, which would not be so uncommon due to the rapid
doubling time of cells?

To offer an insightful response to this intriguing puzzle, we have formulated
a mathematical framework that draws upon analytical components from standard
fitness landscape models, incorporating extensions to account for the effects
of adaptive mutagenesis and ecological influence. Our aim is geared toward
providing a deeper quantitative understanding of stress-induced mutagenesis.
With the model we have developed, not only do we reaffirm the conventional
wisdom that an increase in exploration via a uniformly elevated mutation rate
decreases population success as derived in Eq. (11), but we also show that
focusing on enhancing the mutation rate of an unfit phenotype sub-population –
the essential feature of stress-induced mutagenesis observed in living systems
Cirz (2005); Bos (2015) – can significantly and universally boost up
population success as implied by Eq. (16) and Fig. 2.

Stress-induced mutagenesis drives the evolution of bacterial pathogens,
promoting the emergence of dangerous viral variants and antibiotic resistance
Taddei (1997). This carries significant implications for infectious disease
and biology. As an illustrative example in recent times, stress-induced
mutagenesis has been discovered to have a notable impact on the evolution of
the SARS-CoV-2 virus and the emergence of novel variants Kemp (2021).
Furthermore, antibiotic misuse has led to the development of antibiotic-
resistant bacteria, and this evolutionary mechanism has accelerated the
process by promoting mutations that confer resistance Cirz (2005).
Additionally, stress-induced mutagenesis may contribute to cancer treatment
resistance, highlighting the need for further studies Fitzgerald (2017). Our
paper has provided a mathematical framework for comprehending stress-induced
mutagenesis, which, as we have demonstrated, can potentially unveil novel and
unforeseen emergent population dynamics. Including better representation for
ecological influence in the framework here, beside just a mean-field carrying
capacity, is also a good territory to explore Wang (2021); Phan (2021). By
laying down a rigorous quantitative foundation, we can probe further into the
behavior of adaptive mutagenesis under diverse scenarios, proffering rigorous
insights into this curious sector of evolution. We anticipate further
investigations of our approach to deepen the understanding of stress-induced
mutagenesis, such as why a stochastic dynamical environment can effectively
target highly mutating agents and whether that can be utilized for pathogenic
and neoplastic treatments Day (2022); Wang (2022); Levine (2022).

##  Appendix A The Estimation of $\omega_{l}/\mathscrsfs{D}$ as a Function of
$\epsilon$

For each natural dimensionality $\mathscrsfs{D}=1,2,3,...$, to find the value
$\omega_{l}$ for every given $\epsilon$, we numerically investigate the LHS
and RHS of Eq. (15) as functions of $\omega_{l}$ and estimate their
intersection inside the range $[0,\mathscrsfs{D}]$.

The code for this estimation is available from the corresponding author upon
request. We have executed it on MatLab 2020b.

##  Appendix B Some Calculations for Approximations and Asymptotes of
$\omega_{l}\left(\mathcal{D},\epsilon\right)$

As mentioned, we are interested in finding the solution of $\omega_{l}$ for
$\epsilon\gg 1$, since this is the parameter regime observed in the stress-
induced mutagenesis of real living bacteria Cirz (2005); Bos (2015). By noting
that $\omega_{h}\to 0$ as $\epsilon\to\infty$, we use the following confluent
hypergeometric expansion:

$\operatorname{U}\left(a,b,q\right)\xrightarrow{q\rightarrow
0}\frac{\Gamma\left(1-b\right)}{\Gamma\left(a-b+1\right)}+\frac{\Gamma\left(b-1\right)}{\Gamma\left(a\right)}q^{1-b},$
(17)

to approximate the RHS of Eq. (15). This expansion exhibits different
behaviors in three different cases: $\mathscrsfs{D}=1$, $\mathscrsfs{D}=2$,
$\mathscrsfs{D}\geq 3$.

For $\mathscrsfs{D}=1$ and $\mathscrsfs{D}=2$, we see from Appendix A that the
value of $\omega_{l}$ gradually approaches 0 as $\epsilon\to\infty$. For this
reason we would like to expand the LHS of Eq. (15) in the limit $\omega_{l}\to
0^{+}$. The leading term in this expansion is:

$2\alpha(\mathscrsfs{D})\omega_{l}^{2}$ (18)

where:

$\begin{split}\alpha(\mathscrsfs{D})=&\frac{\frac{d}{d\zeta}L\left(\zeta,\frac{\mathscrsfs{D}}{2},0\right)\Big{|}_{\zeta=-\frac{\mathscrsfs{D}}{4}-1}-4\operatorname{L}\left(-\frac{\mathscrsfs{D}}{4}-2,\frac{\mathscrsfs{D}}{2}+1,0\right)}{4\operatorname{L}\left(-\frac{\mathscrsfs{D}}{4},\frac{\mathscrsfs{D}}{2}-1,0\right)}\\\
&\ \ \ \ \ \ \ \
-\dfrac{\operatorname{L}\left(-\frac{\mathscrsfs{D}}{4}-1,\frac{\mathscrsfs{D}}{2},0\right)\left[\frac{d}{d\zeta}L\left(\zeta,\frac{\mathscrsfs{D}}{2}-1,0\right)\Big{|}_{\zeta=-\frac{\mathscrsfs{D}}{4}}-4\operatorname{L}\left(-\frac{\mathscrsfs{D}}{4}-1,\frac{\mathscrsfs{D}}{2},0\right)\right]}{4\operatorname{L}^{2}\left(-\frac{\mathscrsfs{D}}{4},\frac{\mathscrsfs{D}}{2}-1,0\right)}\
.\end{split}$ (19)

For large $\mathscrsfs{D}$ ($\mathscrsfs{D}\geq 3$), the value of $\omega_{l}$
converges to finite constants near $\mathscrsfs{D}^{-}$, which means we can
take asymptotic expansion for $(\omega_{l}-\mathscrsfs{D})$. For that reason,
the LHS of Eq. (15) becomes:

$\omega_{l}\left[1+\frac{1}{2}(\omega_{l}-\mathscrsfs{D})\frac{d}{d\zeta}\operatorname{L}\left(\zeta,\frac{\mathscrsfs{D}}{2},\mathscrsfs{D}\right)\Big{|}_{\zeta=-1}\right].$
(20)

###  B.1 Estimation of $\omega_{l}(\epsilon)$ for $\mathscrsfs{D}=1$

For $\mathscrsfs{D}=1$, following Eq. (21) we have:

$\frac{\operatorname{U}\left(\frac{\mathscrsfs{D}-\omega_{h}}{4}+1,\frac{\mathscrsfs{D}}{2}+1,\omega_{h}\right)}{\operatorname{U}\left(\frac{\mathscrsfs{D}-\omega_{h}}{4},\frac{\mathscrsfs{D}}{2},\omega_{h}\right)}\Bigg{|}_{\mathscrsfs{D}=1}\xrightarrow{\omega_{h}\to
0}\frac{\Gamma\left(3/4\right)}{\Gamma\left(5/4\right)}\omega_{h}^{-1/2}\ .$
(21)

Together with Eq. (18), Eq. (15) thus becomes:

$2\alpha(1)\omega_{l}^{2}\approx\frac{\mathscrsfs{D}}{2}\left[\frac{\Gamma\left(3/4\right)}{\Gamma\left(5/4\right)}\omega_{h}^{-1/2}\right]\omega_{h}\
\ \Longrightarrow\ \
\omega_{l}(\epsilon)=\left[\frac{1}{4\alpha(1)}\frac{\Gamma\left(3/4\right)}{\Gamma\left(5/4\right)}\right]^{2/3}\epsilon^{-1/6}\propto\epsilon^{-1/6}\
.$ (22)

###  B.2 Estimation of $\omega_{l}(\epsilon)$ for $\mathscrsfs{D}=2$

For $\mathscrsfs{D}=2$, we use Eq. (21) to obtain:

$\frac{\operatorname{U}\left(\frac{\mathscrsfs{D}-\omega_{h}}{4}+1,\frac{\mathscrsfs{D}}{2}+1,\omega_{h}\right)}{\operatorname{U}\left(\frac{\mathscrsfs{D}-\omega_{h}}{4},\frac{\mathscrsfs{D}}{2},\omega_{h}\right)}\Bigg{|}_{\mathscrsfs{D}=2}\xrightarrow{\omega_{h}\to
0}-\frac{2}{\omega_{h}\ln\omega_{h}}\ .$ (23)

Plugging this in the LHS and Eq. (18) in the RHS of Eq. (15), we can arrive
at:

$2\alpha(2)\omega_{l}^{2}=\omega_{h}\left(1-\frac{\mathscrsfs{D}}{\omega_{h}\ln\omega_{h}}\right)\
\ \Longrightarrow\ \
\frac{\omega_{l}^{2}}{\epsilon}\ln\frac{\omega_{l}^{2}}{\epsilon}=-\frac{2}{\alpha(2)\epsilon}$
(24)

The solution of this equation can be represented by the Lambert-W function,
which application has been found in protein physics Levine (2021):

$\omega_{l}=\sqrt{\epsilon\exp\left[W_{-1}\left(-\frac{2}{\alpha(2)\epsilon}\right)\right]}\xrightarrow{\epsilon\rightarrow
0}0$ (25)

###  B.3 Estimation of $\omega_{l}(\epsilon)$ for $\mathscrsfs{D}\geq 3$

For $\mathscrsfs{D}\geq 3$, following Eq. (21) we arrive at:

$\frac{\operatorname{U}\left(\frac{\mathscrsfs{D}-\omega_{h}}{4}+1,\frac{\mathscrsfs{D}}{2}+1,\omega_{h}\right)}{\operatorname{U}\left(\frac{\mathscrsfs{D}-\omega_{h}}{4},\frac{\mathscrsfs{D}}{2},\omega_{h}\right)}\Bigg{|}_{\mathscrsfs{D}\geq
3}\xrightarrow{\omega_{h}\to
0}\frac{2(\mathscrsfs{D}-2)}{\mathscrsfs{D}}\omega_{h}^{-1}\ .$ (26)

This and Eq. (20) simplify Eq. (15):

$\omega_{l}\left[1+\frac{1}{2}(\omega_{l}-\mathscrsfs{D})\frac{d}{d\zeta}\operatorname{L}_{\zeta}^{\mathscrsfs{D}/2}\left(\mathscrsfs{D}\right)\Big{|}_{\zeta=-1}\right]=\mathscrsfs{D}-2\
.$ (27)

This equation is quadratic in $\omega_{l}$, which can be solved to obtain:

$\omega_{l}=\left(\frac{\mathscrsfs{D}}{2}-\frac{1}{Y}\right)\left[1+\sqrt{1+\frac{2(\mathscrsfs{D}-2)Y}{\left(1-\frac{\mathscrsfs{D}}{2}Y\right)^{2}}}\right]\
,\
Y=\frac{d}{d\zeta}\operatorname{L}\left(\zeta,\mathscrsfs{D}/2,\mathscrsfs{D}\right)\Big{|}_{\zeta=-1}\
.$ (28)

This gives the asymptotic limit $\omega_{l}\to\mathscrsfs{D}$ as the number of
traits becomes extremely large $\mathscrsfs{D}\to\infty$.

These analysis has been investigated extensively on Mathematica 13.2.

##  Appendix C The Simulation of Population Distribution on the Landscape

We employ a random-walk agent-based simulation of stress-induced mutagenesis
evolution on a $\mathscrsfs{D}=1$ dimensional landscape, in which
$x\in[-L,+L]$ and the large-scale cut-off is at $L=20$. We investigate the
maximum growth rate $R_{0}=1$, the fit region of size $\lambda=1$, the lower
diffusivity $D_{l}=1/6$. The time resolution of our simulation is $dt=0.01$
and we ran a total time $T=100$, where the stationary state has always been
reached by the end of the first half of the run. The simulation data are
obtained by averaging the recorded every time interval $\Delta t=1$ at the
later half of the run.

The number of agents are dynamical, an agent can be added or removed randomly
with some probability $|p(x,t)|$ per dt, adding when $p(x,t)$ is positive,
removing when $p(x,t)$ is negative. This $p(x,t)$ value is associated with the
growth rate $G(x,t)$ at its position $x$ on the landscape. The carrying
capacity is $K=10^{5}$ and we start with $B(t=0)=10^{4}$ agents randomly
distributed:

$p(x,t)=G(x,t)dt\ ,\
G(x,t)=R_{0}\left[1-\left(\frac{x}{\lambda}\right)^{2}\right]\left[1-\frac{B(t)}{K}\right]\
.$ (29)

The walking step per $dt$ of an agent is selected randomly from a Gaussian
distribution centered at $0$ with a standard deviation $W(x,t)$ associated
with the diffusivity $D[G(x,t)]$ at its position $x$ on the landscape:

$W(x,t)=\sqrt{2D[G(x,t)]dt}\ ,\ D[G]=D_{l}\left[1+\epsilon\Theta(-G)\right]\
.$ (30)

Here we investigate the stress-induced strength $\epsilon=0,3,10$.

The code for this simulation is available from the corresponding author upon
request. We have executed it on MatLab 2020b.

## References

  * Gupta (2006) Gupta, A. K.; Smith, K. G.; and Shalley, C. E.  The interplay between exploration and exploitation.  Academy of management journal 2006, 49(4), 693-706.
  * Uotila (2009) Uotila, J.; Maula, M.; Keil, T.; and Zahra, S. A.  Exploration, exploitation, and financial performance: Analysis of S$\&$P 500 corporations.  Strategic management journal 2009, 30(2), 221-231.
  * Berger-Tal (2014) Berger-Tal, O.; Nathan, J.; Meron, E.; and Saltz, D.  The exploration-exploitation dilemma: a multidisciplinary framework.  PloS one 2014, 9(4), e95693.
  * Guisado-González (2017) Guisado-González, M.; González-Blanco, J.; and Coca-Pérez, J. L.  Analyzing the relationship between exploration, exploitation and organizational innovation.  Journal of Knowledge Management 2017.
  * Alba (2005) Alba, E.; and Dorronsoro, B.  The exploration/exploitation tradeoff in dynamic cellular genetic algorithms.  IEEE transactions on evolutionary computation 2005, 9(2), 126-142.
  * Wright (1932) Wright, S.  The roles of mutation, inbreeding, crossbreeding, and selection in evolution 1932, 355-366.
  * Kimura (1962) Kimura, M.  On the probability of fixation of mutant genes in a population.  Genetics 1962, 47(6), 713.
  * Kimura (1964) Kimura, M.  Diffusion models in population genetics.  Journal of Applied Probability 1964, 1(2), 177-232.
  * Phan (2021) Phan, T. V.; Wang, G.; Do, T. K.; Kevrekidis, I. G.; Amend, S.; Hammarlund, E.; Pienta, K.; Brown, J; Liu, L.; and Austin, R. H.  It doesn’t always pay to be fit: Success landscapes.  Journal of Biological Physics 2021, 47, 387-400.
  * Phan (2020) Phan, T. V.; Morris, R.; Black, M. E.; Do, T. K.; Lin, K. C.; Nagy, K.; Sturm, J.C.; Bos, J.; and Austin, R. H.  Bacterial route finding and collective escape in mazes and fractals.  Physical Review X 2020, 10(3), 031017.
  * Bjedov (2003) Bjedov, I.; Tenaillon, O.; Gerard, B.; Souza, V.; Denamur, E.; Radman, M.; Taddei, F.; and Matic, I.  Stress-induced mutagenesis in bacteria.  Science 2003, 300(5624), 1404-1409.
  * Fitzgerald (2017) Fitzgerald, D. M.; Hastings, P. J.; and Rosenberg, S. M.  Stress-induced mutagenesis: implications in cancer and drug resistance.  Annual Review of Cancer Biology 2017, 1, 119-140.
  * Imasheva (1999) Imasheva, A.G.  Environmental stress and genetic variation in animal populations.  Genetika 1999, 35, 421–431.
  * Hoffmann (2000) Hoffmann, A.A.; and Hercus, M.J.  Environmental stress as an evolutionary force.  Bioscience 2000, 50, 217–226.
  * Cirz (2005) Cirz, R. T.; Chin, J. K.; Andes, D. R.; de Crécy-Lagard, V.; Craig, W. A.; and Romesberg, F. E.  Inhibition of mutation and combating the evolution of antibiotic resistance.  PLoS biology 2005, 3(6), e176.
  * Bos (2015) Bos, J.; Zhang, Q.; Vyawahare, S.; Rogers, E.; Rosenberg, S. M.; and Austin, R. H.  Emergence of antibiotic resistance from multinucleated bacterial filaments.  Proceedings of the National Academy of Sciences 2015, 112(1), 178-183.
  * Foster (2007) Foster, P.L.  Stress-induced mutagenesis in bacteria.  Critical reviews in biochemistry and molecular biology 2007, 42(5), 373-397.
  * Pribis (2022) Pribis, J. P.; Zhai, Y., Hastings; P. J.; and Rosenberg, S. M.  Stress-Induced Mutagenesis, Gambler Cells, and Stealth Targeting Antibiotic-Induced Evolution.  Mbio 2022, 13(3), e01074-22.
  * Zhang (2011) Zhang, Q.; Lambert, G.; Liao, D.; Kim, H.; Robin, K.; Tung, C. K.; Pourmand, N.; and Austin, R. H.  Acceleration of emergence of bacterial antibiotic resistance in connected microenvironments.  Science 2011, 333(6050), 1764-1767.
  * Wu (2014) Wu, A.; Liao, D.; Tlsty, T. D.; Sturm, J. C.; and Austin, R. H.  Game theory in the death galaxy: interaction of cancer and stromal cells in tumour microenvironment.  Interface Focus 2014, 4(4), 20140028.
  * Li (2021) Li, J.; Phulpoto, I. A.; Zhang, G.; and Yu, Z.  Acceleration of emergence of E. coli antibiotic resistance in a simulated sublethal concentration of copper and tetracycline co-contaminated environment.  AMB Express 2021, 11, 1-11.
  * Murray (2022) Murray, C.J.L.; Ikuta, K.S.; Sharara, F.; Swetschinski, L.; Robles-Aguilar, G.; Gray, A.; Han, C.;Bisignano, C.; Rao, P.; Wool, E.; Johnson, S.C.; Browne, A.J.; Chipeta, M.G.; Fell, F.; Hackett, S.; Haines-Woodhouse, G.; Kashef-Hamadani, B.H.; Kumaran, E.A.P.; McManigal B.; Agarwal, R.; Akech, S.; Albertson, S.; Amuasi, J.; Andrews, J.; Aravkin, A.; Ashley, E.; Bailey, F.; Baker, S.; Basnyat, B.; Bekker, A.; Bender, R.; Bethou A.; Bielicki, J.; Boonkasidecha, S.; Bukosia, J.; Carvalheiro, C.; Castañeda-Orjuela, C.; Chansamouth, V.; Chaurasia, S.; Chiurchiù, S; Chowdhury, F.; Cook, A.J.; Cooper, B.; Cressey, T.R.; Criollo-Mora, E.; Cunningham, M.; Darboe, S.; Day, N.P.J.; De-Luca, M.; Dokova, K.; et al.  Global burden of bacterial antimicrobial resistance in 2019: a systematic analysis.  Lancet 2022, 399, 629–655.
  * Ram (2014) Ram, Y.; and Hadany, L.  Stress-induced mutagenesis and complex adaptation.  Proceedings of the Royal Society B: Biological Sciences 2014, 281(1792), 20141025.
  * Tsoularis (2002) Tsoularis, A.; and Wallace, J.  Analysis of logistic growth models.  Mathematical biosciences 2002, 179(1), 21-55.
  * Gavrilets (2018) Gavrilets, S.  Fitness landscapes and the origin of species (MPB-41); Publisher: Princeton University Press, 2018.
  * Vincent (2005) Vincent, T. L; and Brown, J. S.  Evolutionary game theory, natural selection, and Darwinian dynamics.; Publisher: Cambridge University Press, 2005.
  * Phan (2018) Phan, T. V.; Morris, R. J.; Lam, H. T.; Hulamm, P.; Black, M. E.; Bos, J.; and Austin, R. H.  Emergence of Escherichia coli critically buckled motile helices under stress.  Proceedings of the National Academy of Sciences 2018, 115(51), 12979-12984.
  * Risken (1996) Risken, H.; and Caugheyz, T.K.  Fokker-planck equation: Methods of Solution and Applications; Publisher: Springer Berlin Heidelberg 1996.
  * Getz (1991) Getz, W.M.  A unified approach to multispecies modeling.  Natural Resource Modeling 1991, 5(4), 393-421.
  * Getz (1994) Getz, W.M.  A metaphysiological approach to modeling ecological populations and communities.  in Frontiers in Mathematical Biology, 411-442; Publisher: Springer Berlin Heidelberg, 1994.
  * Richards (1959) Richards, F.J.  A flexible growth function for empirical use.  Journal of experimental Botany 1959, 10(2), 290-301.
  * Swartz (2022) Swartz, D. W.; Ottino-Löffler, B.; and Kardar, M.  Seascape origin of Richards growth.  Physical Review E 2022, 105(1), 014417.
  * Tran (2022) Tran, T.; and Kardar, M.  Extinction and growth on an inhomogeneous Seascape.  Turkish Journal of Physics 2022, 46(6), 173-187.
  * McCandlish (2011) McCandlish, D. M.  Visualizing fitness landscapes.  Evolution 2011, 65(6), 1544-1558.
  * Heaviside (1893) Heaviside, O.  On the electromagnetic effects due to the motion of electrification through a dielectric.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 1893, 183(0), 423-480.
  * Abramowitz (1948) Abramowitz, M.; and Stegun, I. A.  FHandbook of mathematical functions with formulas, graphs, and mathematical tables (Vol. 55); Publisher: US Government printing office.
  * Wang (2022) Wang, G.; Phan, T. V.; Li, S.; Wang, J.; Peng, Y.; Chen, G.; Goldman, D.I.; Levin, S.A.; Austin, R.H.; Liu, L.  Robots as models of evolving systems.  Proceedings of the National Academy of Sciences 2022, 119(12), e2120019119.
  * Dobzhansky (1973) Dobzhansky, T.  Nothing in biology makes sense except in light of evolution.  Am.Biol.Teach. 1973, 35, 125-9.
  * Morris (2017) Morris, R. J.; Phan, T. V.; Black, M.; Lin, K. C.; Kevrekidis, I. G.; Bos, J. A.; and Austin, R. H.  Bacterial population solitary waves can defeat rings of funnels.  New Journal of Physics 2017, 19(3), 035002.
  * NRC (1990) Board, Space Studies, and National Research Council.  The Search for Life’s Origins: Progress and Future Directions in Planetary Biology and Chemical Evolution; Publisher: National Academies Press, 1990.
  * Taddei (1997) Taddei, F.; Radman, M.; Maynard-Smith, J.; Toupance, B.; Gouyon, P. H.; and Godelle, B.  Role of mutator alleles in adaptive evolution.  Nature 1997, 387(6634), 700-702.
  * Kemp (2021) Kemp, S.A.; Collier, D.A.; Datir, R.P.; Ferreira, I.A.; Gayed, S.; Jahun, A.; et al  SARS-CoV-2 evolution during treatment of chronic infection.  Nature 2021, 592(7853), 277-282.
  * Wang (2021) Wang, G.; Phan, T. V.; Li, S.; Wombacher, M.; Qu, J.; Peng, Y.; Chen, G.; Goldman, D.I.; Levin, S.A.; Austin, R.H.; and Liu, L.  Emergent field-driven robot swarm states.  Physical review letter 2021, 126(10), 108002.
  * Phan (2021) Phan, T. V.; Wang, G.; Liu, L.; and Austin, R. H.  Bootstrapped motion of an agent on an adaptive resource landscape.  Symmetry 2021, 13(2), 225.
  * Day (2022) Day, C.  Evolving robots could optimize chemotherapy.  Physics Today 2022, 75(3).
  * Levine (2022) Levine, H.  Let the robotic games begin.  Proceedings of the National Academy of Sciences 2022, 119(17), e2204152119.
  * Levine (2021) Khuri, R.R.; Phan, T.V.; and Austin, R.H.  Protein dynamics implications of the low-and high-temperature denaturation of myoglobin.  Physical Review E 2021, 104(3), 034414.


# Learning to Sail Dynamic Networks:
The MARLIN Reinforcement Learning Framework for Congestion Control in Tactical
Environments

Raffaele Galliera12, Mattia Zaccarini3, Alessandro Morelli1, Roberto
Fronteddu1,
Filippo Poltronieri3, Niranjan Suri14, Mauro Tortonesi3  1 Florida Institute
for Human & Machine Cognition (IHMC)
Email: {rgalliera, amorelli, rfronteddu<EMAIL_ADDRESS>2 Department of
Intelligent Systems & Robotics - The University of West Florida (UWF),
Pensacola, FL, USA
Email<EMAIL_ADDRESS>3 Distributed Systems Research Group,
University of Ferrara, Ferrara, Italy
Email: {mattia.zaccarini, filippo.poltronieri<EMAIL_ADDRESS>4 US
Army DEVCOM Army Research Laboratory (ARL), Adelphi, MD, USA
Email<EMAIL_ADDRESS>
###### Abstract

Conventional Congestion Control (CC) algorithms, such as TCP Cubic, struggle
in tactical environments as they misinterpret packet loss and fluctuating
network performance as congestion symptoms. Recent efforts, including our own
MARLIN, have explored the use of Reinforcement Learning (RL) for CC, but they
often fall short of generalization, particularly in competitive, unstable, and
unforeseen scenarios. To address these challenges, this paper proposes an RL
framework that leverages an accurate and parallelizable emulation environment
to reenact the conditions of a tactical network. We also introduce refined RL
formulation and performance evaluation methods tailored for agents operating
in such intricate scenarios. We evaluate our RL learning framework by training
a MARLIN agent in conditions replicating a bottleneck link transition between
a Satellite Communication (SATCOM) and an UHF Wide Band (UHF) radio link.
Finally, we compared its performance in file transfer tasks against
Transmission Control Protocol (TCP) Cubic and the default strategy implemented
in the Mockets tactical communication middleware. The results demonstrate that
the MARLIN RL agent outperforms both TCP and Mockets under different
perspectives and highlight the effectiveness of specialized RL solutions in
optimizing CC for tactical network environments.

###### Index Terms:

Machine Learning, Congestion Control, Reinforcement Learning, Emulated
Networks, Tactical Networks.

##  I Introduction

Tactical networks present unique challenges for CC. While operating in
demanding circumstances marked by limited resources, unreliable links,
frequent disconnections, and varying levels of connectivity, these networks
need to support critical real-time functionalities to facilitate mission
applications such as Command, Control, Communications, Computers, Cyber,
Intelligence, Surveillance, and Reconnaissance (C5ISR).

Traditional CC algorithms, such as those embedded in transport protocols like
TCP, face significant difficulties in maintaining efficient communications
within tactical networks, as they were developed for wired environments and
misinterpret as congestion symptoms phenomena such as packet losses and
temporary unreachability, that are very common in tactical networks, thus
severely and unnecessarily reducing transmission speeds. Optimizing CC in
unreliable networks necessitates innovative approaches that are able to cope
with their dynamic and unpredictable nature.

In this context, the combination of RL techniques with Deep Learning (DL) for
policy parameterization, often referred to as Deep Reinforcement Learning, has
emerged as a promising approach. Deep RL demonstrated remarkable robustness
across diverse domains, and its application in computer networks offers new
possibilities for addressing the challenges of CC [1, 2]. By leveraging RL,
agents can be trained to learn optimal policies through interactions with the
network environment, enabling more efficient and reliable communications.

Despite numerous efforts showing promising results, the performance of these
RL agents still falls short of generalization capabilities, especially when
unreliable and unpredictable links are encountered. This performance gap can
be attributed to various factors either linked to the learning problem itself,
such as designing the environment where the agent resides and its "learning
curriculum", or to the challenging and partially observable nature of the
networking environment. Addressing these challenges requires the development
of solutions that enable RL agents to learn effective policies for CC
decision-making while experiencing complex and dynamic scenarios. However,
creating dedicated environments for conducting experiments in such scenarios
can be challenging, costly, and even infeasible in real network environments.
To overcome this last constraint, the importance of an accurate emulation
environment cannot be overstated.

This paper extends our work within the MARLIN project [3] by proposing an RL
framework that leverages an accurate and parallelizable emulation environment
to reenact the conditions of a tactical network, thus allowing agents to
experience a wide range of dynamic behaviors. To better evaluate the decision-
making effectiveness of our agent, we also introduce a novel metric, the RTT
Transition Impact (RTI), based on the maximum Round-Trip Time (RTT) detected
during a communication involving a link transition, that allows to measure the
agent’s responsiveness to link changes and its queue management capabilities.

To evaluate our learning framework, we trained an RL agent in an emulated
environment replicating a bottleneck link transition between a SATCOM link and
an UHF radio link - rather typical in tactical networks. We then evaluated the
agent performance in a file transfer task and compared it against conventional
CC algorithms such as TCP Cubic [4], as well as CC algorithms implemented
within communication middlewares tailored for tactical environments such as
Mockets [5], which the same MARLIN utilizes as partnering protocol. The
results demonstrate that the exposure of the agent to complex networking
scenarios enables training policies able to achieve competitive decision-
making performance, and validate the RL training approach based on accurate
emulation and purposely designed valuation metrics for specialized
environments like tactical networks.

##  II Congestion Control in Tactical Networks

Tactical networks represent a unique communications environment due to a
combination of particularly harsh and dynamic network conditions paired with a
high degree of heterogeneity in the network technologies employed and,
consequently, the characteristics of the links. This situation demands
efficient data transmission strategies that can quickly respond to changes in
the status of the network, which makes the role of transport protocols and CC
algorithms of crucial importance in tactical networks and especially
challenging.

Traditional and de-facto standard CC algorithms for TCP, such as CUBIC, were
not designed for this type of scenario and, as a consequence, tend to
underperform significantly or, in some cases, cause connections to break
altogether, leading to significant loss of network resources [6]. Frequent
packet loss causes loss-based CC algorithms, which interpret loss as a sign of
congestion, to step out of the slow start phase and decrease their congestion
window, impeding the effectiveness of the bandwidth discovery process. In
contrast, delay-based CC algorithms typically misinterpret abrupt changes in
the end-to-end latency, such as those caused by changes in the routing rules,
as a sign that they are allowed to either increase or reduce their Congestion
Window (CWND) significantly. This behavior can again lead to severe
underutilization of already scarce network resources.

This fostered the design of dedicated solutions to overcome the limitations of
TCP in tactical environments, such as Mockets [5], a User Datagram Protocol
(UDP)-based communication middleware that allows extensive customization of
its configuration parameters. By enabling adjustments to parameters like the
CWND, the Pending Packet Queue Size, and the Selective Acknowledgment
Transmission Timeout, Mockets offers the flexibility to fine-tune its settings
according to the specific requirements of different communication channels.
Mockets implements a very aggressive CC algorithm whose purpose is to fully
utilize the link bandwidth in the presence of variable communication latency
and elevated packet loss. While it presents several advantages over TCP and
other transport protocols in degraded environments [6], the existing CC fails
to share link capacity with other communication flows going through the same
links and cannot adapt quickly to changes in the available bandwidth.

These limitations are especially visible in case of heterogeneous network
technologies. In fact, abrupt changes from high-bandwidth, high-latency links
like SATCOM to a low-bandwidth, low-latency links like UHF, and back, can
wreak havoc on the performance of CC solutions. To compound these issues,
tactical environments exhibit: routing rules that can change during a mission;
highly varying wireless performance; node mobility; and the enemy’s activity,
which could jam communications and/or render nodes unreachable. These
challenging scenarios necessitate the utilization of advanced techniques
capable of enabling dynamic and responsive reconfiguration in response to
evolving network conditions.

##  III Leveraging Container Networks to Emulate Challenging Networking
Scenarios

Figure 1: The MARLIN containerized training framework.

Complex networking scenarios, such as those encountered in tactical
environments, are known for their intricacy, often involving sets of
heterogeneous entities with distinct behaviors [7]. However, RL environments
dedicated to CC are usually not designed to capture such complexities in their
training setting, making the agents prone to underperform when deployed in
these scenarios.

To bridge this gap, our work enhances the MARLIN framework111Code available at
https://github.com/RaffaeleGalliera/marlin-rlcc by incorporating customizable
link behaviors, which in this case we implement to resemble the complexities
of tactical environments. By introducing these diverse scenarios into the
training environment, we can expose the RL agent to a wider range of
experiences, enabling it to learn and adapt to the intricate behaviors
exhibited by the heterogeneous entities in tactical networks.

Specifically, we extend the MARLIN solution to employ a network generator
module. In succinct terms, such a component is responsible to create custom
container networks, define the behavior of the links involved, and generate
background traffic from competitive/collaborative sources. The reader can
refer to Figure 1 for a summary of the framework, which also includes the
topology defined for this work (detailed in V).

The network generator module defines a Remote Procedure Call (RPC) interface
exposing a series of methods bonded to a Containernet [8] instance operating
on the host machine. Containernet is a lightweight container-based network
emulation tool built on Mininet that allows users to create and manage virtual
networks with Docker containers for the purpose of simulating complex network
scenarios.

Once started, our network generator module defines a starting topology running
specific Docker images for every container involved. The RPC methods are then
used by the RL environment to create and define diverse scenarios. Every time
the environment is initialized or reset, a dedicated RPC method is invoked,
setting the links to their original values. In the meantime, the same
environment will start the communication between the Mockets sender and its
receiving counterpart, which effectively starts the collaboration between the
agent and the protocol. Other hosts are prompted to generate traffic,
competing with the RL agent for link usage. These containers will induce
different TCP and UDP traffic behaviors following a Multi-Generator Network
Test Tool (MGEN) [9] script defined through the RPC arguments.

In Section IV-C, we detail how such a controlled setting might guide the agent
toward effective policies.

##  IV A Reinforcement Learning Formulation for Congestion Control in Tactical
Environments

In this section, we present the fundamental RL [10] components on which this
version of MARLIN relies. In particular, we address the challenge of
determining the timing of actions in the tactical scenario and propose a
reward function specifically designed for these complex environments.

###  IV-A Observation and Action space

MARLIN utilizes a state representation consisting of 14 features collected
during the timeframe following the previous action at step $t-1$. The state
space is expanded by computing 7 statistics for each feature, including the
last value, mean, standard deviation, minimum, maximum, Exponential Moving
Average (EMA), and difference from the previous state $s_{t-1}$. The previous
10 states are stacked together, forming an observation history. Observations
are then normalized using a moving average. The action space, instead,
comprises continuous actions within the range of $[-1,1]$, denoting the
percentage gain in congestion window (CWND) size. For example, an action of .3
would increase the CWND by 30%, while -.3 would instead reduce it by the same
quantity.

###  IV-B Timing Decisions

When designing a RL-based CC agent, determining when the next action is going
to be taken constitutes one of the first questions that need to be answered.
In the original work [3], the authors opted for a heuristic based on the
Smoothed Round-Trip Time (SRTT) detected by the partnering protocol. This
approach becomes problematic when the agent is trained within links with
delays in the order of seconds, as it would increase the time needed to
complete the same training process by several orders of magnitude, without any
guarantee of convergence to a policy with the same quality.

To address this limitation, our work utilizes a predefined time window of
100ms. Within this timeframe, the partnering protocol collects the necessary
statistics, which will be used to form the state of the agent.

###  IV-C Rewarding the agent in Tactical Networks

We designed a reward function incorporating penalties based on
retransmissions, proving it to be highly suitable to tactical networks, owing
to a multitude of compelling reasons.

By integrating penalties for retransmissions into the reward function, the RL
agent is actively incentivized to minimize their occurrence, thereby improving
the reliability of data transmission. This behavior curtails potential delays
and data loss, which harmoniously aligns with the specific requirements and
constraints inherent to tactical networks.

Building upon the principles delineated in [3] with the controlled training
scenarios, it is possible to approximate the traffic influx from competing
sources at each time step $t$. Consequently, reward functions can be devised
to both confer greater rewards as the agent progressively approaches the
optimal utilization of the available bandwidth and greater penalties when
retransmissions are present, thereby incentivizing efficient network usage.
These concepts can be merged into a single reward function described as
follows:

$r_{t}=-\frac{target_{t}[1+r*(1-loss_{c})]}{target_{t}+acked_{t}^{cumulative}}$
(1)

where $target_{t}$ is an approximation of the number of bytes the agent could
have delivered up to step $t$ since the beginning of the episode in order to
fully utilize the link and $acked_{t}^{cumulative}$ represents the number of
kilobytes there were acknowledged by the receiver until step $t$. The value
$r$ quantifies the number of retransmissions detected by the partnering
protocol, while $loss_{c}$ is the decimal representation of the current random
packet loss.

It is worth mentioning that a strictly negative reward function promotes the
agent to accumulate the smallest amount of penalties. The penalty received is
smaller the closer the agent gets, at each step, to utilize the link to the
best of its possibilities, avoiding overwhelming the other traffic flows.

((a)) Background traffic during the SATCOM link usage

((b)) Background traffic during the UHF link usage

Figure 2: Characterization of the background traffic patterns.

##  V Experimental Results

Figure 3: Performance comparison in terms of average transfer time and RTI
across the evaluation episodes.

In order to replicate the wide-ranging variations commonly observed in
tactical networks, we have developed a testbed that serves as a training
scenario and implements an instance of the tactical environment described in
Section II.

In this particular scenario, a data rate of 1Mb/s is initially accessible
through a SATCOM link, which exhibits a 500ms delay due to the satellite-
ground station distance. However, as the tactical network operates in a
dynamic environment, it is crucial for it to respond efficiently even in
situations where the SATCOM link becomes unavailable or experiences sudden
performance degradation. To address such behavior, we simulated a transition
to an alternative UHF radio link with a 256Kbps data rate, a reduced delay of
125ms, and an assumed 3% packet loss. Following our framework detailed in
Section III it is possible to reproduce such a dynamic environment by defining
three elements: the network architecture, the link behavior, and the
background traffic patterns.

In Figure 1, we present the design of a dumbbell architecture that separates
two networks, denoted as $LS$ and $RS$. The sender hosts of both Mockets and
the MGEN traffic generator are situated on network $LS$, while their receiving
counterparts reside on network $RS$. It is worth noting that the network
architecture remains unchanged throughout the scenario.

To replicate the desired scenario, we implement the linking behavior within
our network generator module. Initially, the bottleneck link is configured to
emulate the properties of the SATCOM link. However, after a duration of 10
seconds, a link transition is triggered, reconfiguring the link to embody the
characteristics of the UHF radio link. This adjustment accurately reproduces
the intended scenario within the testbed.

In the meantime, background traffic flows are emitted by the competing
sources, following the sequences presented in Figure 2. Each sequence is
repeated every 8 seconds and includes two parallel components divided into
elephant flows, i.e., long-lived data transfers that represent a large
percentage of the total traffic, and mice flows, i.e., short-lived data
transfers at low throughput. Elephant flows are alternated every two seconds,
while two mice flows are continuously generating very short-lived (in the
order of milliseconds) traffic bursts with intervals that follow a Poisson
distribution.

Traffic patterns are adjusted to align with the characteristics of the UHF
link at the time of transition, as illustrated in Figure 2(b). This effect
simulates a perfect adaptation of the background traffic flows to the link
change, assuming that competing traffic sources are able to seamlessly adjust
to link transitions. It is also supposed that the occurrence of significant
UDP elephant traffic, such as video streams, is interrupted from the competing
sources while the communication is routed through the UHF link.

###  V-A Training

When training a RL agent, in addition to its objectives, we also need to
define the “horizon” towards which it is going to optimize its behavior. For
this purpose, many RL problems define terminal states and the rewards
associated with them serve as a signal for successful completion. Once a
terminal state is reached, the current episode ends, and a new episode begins.
In problems like CC it is challenging to define terminal states as, at least
theoretically, a communication instated between two nodes could continue
indefinitely. For such a reason, MARLIN considers CC a time-unlimited task
during training. In these settings, the agent does not optimize its behavior
towards reaching a predefined terminal state while maximizing its return, but
the focus is shifted to maximizing its returns over an indefinite time
horizon.

Partial episodes, each of which was 200 steps long, were used to diversify the
experience of the agent for an approximate duration of 20 seconds per episode,
considering that the interval between actions was set to 100 milliseconds, as
discussed in Section IV-B.

The algorithms and hyperparameters used mostly followed the default values
presented in [3], except for the buffer size and training duration, which were
reduced and set respectively to 250K and 500K steps. The observation space and
the action space are unaffected compared to the previous work, except for the
maximum value that the CWND might assume, here set to 150KB.

###  V-B Evaluation Method

To evaluate the performance of the trained agent, we conducted a series of
file transfers in our experiments. Each evaluation episode involved
transmitting a 600KB payload. The experiment, repeated 400 times for both TCP
and our version of MARLIN, was separated in batches of 100 transfers under
various random packet loss probabilities of the UHF link, ranging from 0% to
3%. It is important to highlight that the agent was exclusively trained using
a random packet loss probability of 3%. The purpose of this evaluation is to
confirm that the agent did not learn any patterns based on the number of
retransmissions induced by the emulated packet loss.

The results were compared together to the performance of a random agent and
Mockets, both performing only on the 3% packet loss scenario, as they do not
adopt any adaptive mechanism. The latter was configured such that the CWND
matched the bandwidth capacity of the link, before and after the transition.
Enabled by our controlled environment, this seamless process eliminates the
necessity of customary protocol reconfiguration, thereby intensifying the
challenge for our agent. Finally, we approximated an ideal fair behavior,
assuming an algorithm able to split the bandwidth available left by non-
collaborative flows (i.e. UDP traffic) with adaptive flows (i.e. TCP flows),
without causing any RTT increase.

Evaluating CC solutions in dynamic and challenging links necessitates careful
consideration in the selection of appropriate metrics. In our evaluation
process, we focused on three key aspects: the file transfer completion time,
the number of retransmissions incurred, and the sensitivity of the agent to
link changes and its queue management capabilities.

While the first two factors can be straightforwardly quantified, we believe
the latter to be more subtle. Simply relying on the average or maximum RTT
measured in the link might seem reasonable at first sight, but they would
either fail to capture buffer loading peaks or lack information about the
specific link where the peak incurred. These observations emphasize the
necessity of employing an able to capture the responsiveness of the agent to
adapt to the new link characteristics.

To address this challenge, we designed a novel metric called the _RTT
Transition Impact (RTI)_. With the objective of assessing the influence of CC
algorithms on buffer loading when transitioning between the links, the RTI
leverages the use of the normalized maximum RTTs measured separately for each
of the $m$ links involved during communication. By considering individual link
measurements, the RTI provides insights into the impact of link changes on the
performance of the CC algorithm and it can be defined as follows:

$RTI=\ln(\frac{\sum^{m}_{i=1}\frac{rtt_{i,max}}{rtt_{i,nom}}}{m})$ (2)

where $m$ identifies the different links involved. The maximum RTT detected
during the timeframe related to link $i$ is denoted as $rtt_{i,max}$, while
the nominal RTT for the same timeframe is denoted as $rtt_{i,nom}$. By
employing the natural logarithm, we can emphasize the magnitude of the impact
for smaller values, while denoting with a value of 0 a flawless performance.

Through the utilization of the RTI metric, we enhance our evaluation
methodology by capturing the intricate dynamics of tactical network
environments and their effect on CC algorithms. This metric contributes to a
more comprehensive understanding of its ability to adapt to evolving link
conditions, thereby facilitating a more accurate assessment of CC performance
in dynamic and challenging scenarios.

###  V-C Results

Figure 3 presents a comparison between MARLIN, Mockets, TCP Cubic, a random
agent, and the ideal fair behavior in terms of the average transfer time and
the average RTI achieved. The results show that our trained agent exhibits a
moderate improvement in file transfer completion time compared to TCP, with an
average performance of 19.3s compared to 22.20s. However, it is noticeably
slower than Mockets, which achieves an average completion time of 10.59s. This
outcome is expected as Mockets is purposefully designed for tactical
environments and adopts an aggressive policy that introduces high overhead on
the network, as we will further discuss.

In terms of average RTI, MARLIN consistently outperforms its competitors
across different scenarios, regardless of varying packet losses and with
minimal RTI fluctuations between the 4 scenarios. On the other hand, TCP Cubic
exhibits a peculiar and variable behavior in terms of RTI performance. It is
observed that lower levels of packet loss result in higher RTI for TCP Cubic,
indicating potential limitations in its queue management capabilities. This
contrast highlights the robustness and stability of MARLIN in maintaining low
RTI values in our scenario, showcasing promising performance for tactical
network environments.

To conclude and support our results, we assessed the average number of
retransmissions per episode involving MARLIN, Mockets, the random agent, and
TCP Cubic when a 3% packet loss was present on the UHF link. Figure 4
illustrates how our agent is able to limit the average retransmissions to
12.22 against the 53.1 for TCP Cubic and 141.93 for Mockets.

Figure 4: The average number of retransmissions during evaluation with a 3%
random packet loss set on the UHF link.

##  VI Related Work

Considerable emphasis has consistently been placed on the optimization of CC
algorithms. This has led to the proposal of numerous heuristics, including
notable examples such as TCP Vegas [11], BBR [12], and TCP CUBIC [4]. However,
most algorithms achieve satisfactory performance in specific scenarios often
demonstrating poor generalization, especially in dynamic and unreliable
networks [13].

In the ever-evolving field of networking, it is not surprising that recent
efforts have focused also on RL-based optimization techniques, aiming to
develop solutions able to adapt to the complexity and heterogeneity of
networking environments. Approaches such as DRL-CC [2] have focused on
training policies in simulated environments, achieving competitive
performance, but demonstrating several challenges when applied to real
networks. Others, such as MVFST-RL [1] focused on training non-blocking CC RL
agents on the application level, overcoming the limitation of non-blocking
solutions. MARLIN [3], followed a similar non-blocking approach, with a
particular emphasis on training and evaluating the agent end-to-end on a real
network.

To the best of our knowledge, this paper presents a two-fold novel
contribution: it is the first RL environment for CC with a focus on tactical
networks and the first exploiting the flexibility of emulated networks with
containerized applications to train and evaluate a RL agent for CC.

##  VII Conclusion and Future Work

This paper presented a step toward comprehensive RL environments capable of
training policies on emulated networks with dynamic link behavior, paving the
way for studying the application of RL to complex networking scenarios, which
might be impractical or costly in real network environments.

Experimental results demonstrated how MARLIN, with a relatively exiguous
training budget, was able to learn an effective CC policy in the presented
scenario, addressing the shortcomings of both the Mockets’ aggressive strategy
and TCP Cubic. Among the other criteria used for evaluation, we have also
introduced the RTI, a metric useful for capturing the sensitivity of the agent
to link changes.

It is crucial to mention that the framework and the results presented here,
albeit we believe them to be promising and compelling, are still limited and
focus on a specific training scenario. A variety of approaches and methods
still need to be explored, implemented, and evaluated within MARLIN, which
will be part of our future work. These include but are not limited to,
parallel environments, more heterogeneous dynamic network setups, self-timed
decisions, and competitive and/or cooperative Multi-Agent Reinforcement
Learning (MARL) settings.

## References

  * [1] V. Sivakumar, O. Delalleau, T. Rocktäschel, A. H. Miller, H. Küttler, N. Nardelli, M. Rabbat, J. Pineau, and S. Riedel, “Mvfst-rl: An asynchronous rl framework for congestion control with delayed actions,” 10 2019. [Online]. Available: http://arxiv.org/abs/1910.04054
  * [2] Z. Xu, J. Tang, C. Yin, Y. Wang, and G. Xue, “Experience-driven congestion control: When multi-path tcp meets deep reinforcement learning,” _IEEE Journal on Selected Areas in Communications_ , vol. 37, no. 6, pp. 1325–1336, 2019\.
  * [3] R. Galliera, A. Morelli, R. Fronteddu, and N. Suri, “Marlin: Soft actor-critic based reinforcement learning for congestion control in real networks,” in _NOMS 2023-2023 IEEE/IFIP Network Operations and Management Symposium_ , 2023, pp. 1–10.
  * [4] S. Ha, I. Rhee, and L. Xu, “Cubic: A new tcp-friendly high-speed tcp variant,” _SIGOPS Oper. Syst. Rev._ , vol. 42, no. 5, p. 64–74, jul 2008\. [Online]. Available: https://doi.org/10.1145/1400097.1400105
  * [5] E. Benvegnù, N. Suri, M. Tortonesi, and T. Esterrich, “Seamless network migration using the mockets communications middleware,” in _2010 - MILCOM 2010 MILITARY COMMUNICATIONS CONFERENCE_ , 2010, pp. 2298–2303.
  * [6] A. Morelli, M. Provosty, R. Fronteddu, and N. Suri, “Performance evaluation of transport protocols in tactical network environments,” in _MILCOM 2019 \- 2019 IEEE Military Communications Conference (MILCOM)_ , 2019, pp. 30–36.
  * [7] N. Suri, A. Hansson, J. Nilsson, P. Lubkowski, K. Marcus, M. Hauge, K. Lee, B. Buchin, L. Mısırhoğlu, and M. Peuhkuri, “A realistic military scenario and emulation environment for experimenting with tactical communications and heterogeneous networks,” in _2016 International Conference on Military Communications and Information Systems (ICMCIS)_ , 2016, pp. 1–8.
  * [8] M. Peuster, H. Karl, and S. van Rossem, “Medicine: Rapid prototyping of production-ready network services in multi-pop environments,” in _2016 IEEE Conference on Network Function Virtualization and Software Defined Networks (NFV-SDN)_ , Nov 2016, pp. 148–153.
  * [9] N. R. L. N. P. E. A. N. P. R. Group, “Multi-generator (mgen) network test tool,” https://www.nrl.navy.mil/Our-Work/Areas-of-Research/Information-Technology/NCS/MGEN/, 2021\.
  * [10] R. S. Sutton and A. G. Barto, _Reinforcement Learning: An Introduction_. Cambridge, MA, USA: A Bradford Book, 2018.
  * [11] L. S. Brakmo, S. W. O’Malley, and L. L. Peterson, “Tcp vegas: New techniques for congestion detection and avoidance,” in _SIGCOMM_ , 1994.
  * [12] N. Cardwell, Y. Cheng, C. S. Gunn, S. H. Yeganeh, and V. Jacobson, “Bbr: Congestion-based congestion control,” _ACM Queue_ , vol. 14, September-October, pp. 20 – 53, 2016. [Online]. Available: http://queue.acm.org/detail.cfm?id=3022184
  * [13] P. Chaudhary and S. Kumar, “Comparative study of tcp variants for congestion control in wireless network,” in _2017 International Conference on Computing, Communication and Automation (ICCCA)_ , 2017, pp. 641–646.


# Sub-barrier capture reactions with 16,18O and 40,48Ca beams

V.V.Sargsyan1, G.G.Adamian1, N.V.Antonenko1, W. Scheid2, and H.Q.Zhang3 1Joint
Institute for Nuclear Research, 141980 Dubna, Russia
2Institut für Theoretische Physik der Justus–Liebig–Universität, D–35392
Giessen, Germany
3China Institute of Atomic Energy, Post Office Box 275, Beijing 102413, China

###### Abstract

Various sub-barrier capture reactions with beams 16,18O and 40,48Ca are
treated within the quantum diffusion approach. The role of neutron transfer in
these capture reactions is discussed. The quasielastic and capture barrier
distributions are analyzed and compared with the recent experimental data.

###### pacs:

25.70.Jj, 24.10.-i, 24.60.-k
Key words: sub-barrier capture, neutron transfer, quantum diffusion approach

##  I Introduction

From the present experimental data the role of the neutron transfer channel in
the capture (fusion) process cannot be unambiguously inferred Bertulani ; Jia
; Kol ; EPJSub ; EPJSub1 . The fusion excitation functions have been recently
measured for the reactions 16O+76Ge and 18O+74Ge at energies near and below
the Coulomb barrier and the fusion barrier distributions have been extracted
from the corresponding excitation functions Jia . The fusion enhancement due
to the positive $Q_{2n}$-value two neutron ($2n$) transfer channel for
18O+74Ge has not been revealed as compared with the reference system 16O+76Ge
Jia . This is very different from the situation for the reactions
40Ca+124,132Sn Kol  and other systems in literature, which show considerable
sub-barrier enhancements. The enhancement appears to be related to the
existence of positive $Q$ values for neutron transfer.

The purpose of this paper is the theoretical explanation of these experimental
observations. Within the quantum diffusion approach EPJSub ; EPJSub1  we try
to answer the question how strong the influence of neutron transfer in sub-
barrier capture (fusion) reactions
18O+74Ge,52,50Cr,94,92Mo,112,114,118,120,124,126Sn and 40,48Ca+124,132Sn. This
study seems to be important for future experiments indicated in Ref. Jia . In
addition, the new structures of the quasielastic and capture barrier
distributions at deep sub-barrier energies will be discussed.

##  II Model

In the quantum diffusion approach EPJSub ; EPJSub1  the collisions of nuclei
are described with a single relevant collective variable: the relative
distance between the colliding nuclei. This approach takes into consideration
the fluctuation and dissipation effects in collisions of heavy ions which
model the coupling with various channels (for example, coupling of the
relative motion with low-lying collective modes such as dynamical quadrupole
and octupole modes of the target and projectile Ayik333 ). We have to mention
that many quantum-mechanical and non-Markovian effects accompanying the
passage through the potential barrier are taken into consideration in our
formalism EPJSub1 . The nuclear deformation effects are taken into account
through the dependence of the nucleus-nucleus potential on the deformations
and mutual orientations of the colliding nuclei. To calculate the nucleus-
nucleus interaction potential $V(R)$, we use the procedure presented in Refs.
EPJSub1 . For the nuclear part of the nucleus-nucleus potential, the double-
folding formalism with the Skyrme-type density-dependent effective nucleon-
nucleon interaction is used. With this approach many heavy-ion capture
reactions at energies above and well below the Coulomb barrier have been
successfully described EPJSub1 . Note that the diffusion models, which include
the quantum statistical effects, were also treated in Refs. Hofman .

Following the hypothesis of Ref. Broglia , we assume that the sub-barrier
capture in the reactions under consideration mainly depends on the possible
two-neutron transfer with the positive $Q_{2n}$-value. Our assumption is that,
just before the projectile is captured by the target-nucleus (just before the
crossing of the Coulomb barrier) which is a slow process, the $2n$-transfer
($Q_{2n}>0$) occurs that can lead to the population of the excited collective
states in the recipient nucleus SSzilner . So, the motion to the $N/Z$
equilibrium starts in the system before the capture because it is
energetically favorable in the dinuclear system in the vicinity of the Coulomb
barrier. For the reactions considered, the average change of mass asymmetry is
related to the two-neutron transfer. In these reactions the $2n$-transfer
channel is more favorable than $1n$-transfer channel ($Q_{2n}>Q_{1n}$). Since
after the $2n$-transfer the mass numbers, the deformation parameters of the
interacting nuclei, and, correspondingly, the height $V_{b}=V(R_{b})$ of the
Coulomb barrier are changed, one can expect an enhancement or suppression of
the capture. This scenario was verified in the description of many reactions
EPJSub1 .

##  III Results of calculations

All calculated results are obtained with the same set of parameters as in Ref.
EPJSub1  and are rather insensitive to the reasonable variation of them
EPJSub1 . Realistic friction coefficient in the momentum $\hbar\lambda$=2 MeV
is used which is close to those calculated within the mean field approaches
obzor . The parameters of the nucleus-nucleus interaction potential $V(R)$ are
adjusted to describe the experimental data at energies above the Coulomb
barrier corresponding to spherical nuclei. The absolute values of the
quadrupole deformation parameters $\beta_{2}$ of even-even deformed nuclei are
taken from Ref. Ram . In Ref. Ram , the quadrupole deformation parameters
$\beta_{2}$ are given for the first excited 2+ states of nuclei. For the
nuclei deformed in the ground state, the $\beta_{2}$ in 2+ state is similar to
the $\beta_{2}$ in the ground state and we use $\beta_{2}$ from Ref. Ram  in
the calculations. For the double magic nucleus 16O, in the ground state we
take $\beta_{2}=0$. Since there are uncertainties in the definition of the
values of $\beta_{2}$ in light- and medium-mass nuclei, one can extract the
quadrupole deformation parameters of these nuclei from a comparison of the
calculated capture cross sections with the existing experimental data. By
describing the reactions 18O+208Pb, where there are no neutron transfer
channels with positive $Q$-values, we extract $\beta_{2}=0.1$ for the ground-
state of 18O EPJSub1 . This extracted value is used in our calculations.

###  III.1 Effect of neutron transfer in reactions with beams 40,48Ca

To eliminate the influence of the nucleus-nucleus potential on the capture
(fusion) cross section and to make conclusions about the role of deformation
of colliding nuclei and the nucleon transfer between interacting nuclei in the
capture (fusion) cross section, a reduction procedure is useful Gomes . It
consists of the following transformations:

$E_{\rm c.m.}\rightarrow x=\dfrac{E_{\rm
c.m.}-V_{b}}{\hbar\omega_{b}},\qquad\sigma_{cap}\rightarrow\sigma_{cap}^{red}=\dfrac{2E_{\rm
c.m.}}{\hbar\omega_{b}R_{b}^{2}}\sigma_{cap},$

where $\sigma_{cap}=\sigma_{cap}(E_{\rm c.m.})$ is the capture cross section
at bombarding energy $E_{\rm c.m.}$. The frequency
$\omega_{b}=\sqrt{V^{{}^{\prime\prime}}(R_{b})/\mu}$ is related with the
second derivative $V^{{}^{\prime\prime}}(R_{b})$ of the total nucleus-nucleus
potential $V(R)$ (the Coulomb + nuclear parts) at the barrier radius $R_{b}$
and the reduced mass parameter $\mu$. With these replacements we compared the
reduced calculated capture (fusion) cross sections $\sigma_{cap}^{red}$ for
the reactions 40,48Ca+124,132Sn (Fig. 1). The choice of the projectile-target
combination is crucial, and for the systems studied one can make unambiguous
statements regarding the neutron transfer process with a positive $Q$-value
when the interacting nuclei are double magic or semi-magic spherical nuclei.
In this case one can disregard the strong direct nuclear deformation effects.
In Fig. 1, one can see that the reduced capture cross sections in the
reactions 40Ca+124,132Sn with the positive $Q_{2n}$-values strongly deviate
from those in the reactions 48Ca+124,132Sn, where the neutron transfers are
suppressed because of the negative $Q$-values.

Figure 1: (Colour online) The calculated reduced capture cross sections versus
$(E_{\rm c.m.}-V_{b})/(\hbar\omega_{b})$ in the reactions 40Ca+124Sn (solid
line), 48Ca+124Sn (dashed line), 48Ca+124Sn (dotted line), and 48Ca+132Sn
(dash-dotted line).

After two-neutron transfer in the reactions
40Ca($\beta_{2}=0$)+124Sn($\beta_{2}=0.1$)$\to^{42}$Ca($\beta_{2}=0.25$)+122Sn($\beta_{2}=0.1$)
($Q_{2n}$=5.4 MeV) and
40Ca($\beta_{2}=0$)+132Sn($\beta_{2}=0$)$\to^{42}$Ca($\beta_{2}=0.25$)+130Sn($\beta_{2}=0$)
($Q_{2n}$=7.3 MeV) the deformation of the light nucleus increases and the mass
asymmetry of the system decreases and, thus, the value of the Coulomb barrier
decreases and the capture cross section becomes larger (Fig. 1). So, because
of the transfer effect the systems 40Ca+124,132Sn show large sub-barrier
enhancements with respect to the systems 48Ca+124,132Sn. We observe that the
$\sigma_{cap}^{red}$ in the 40Ca+124Sn (48Ca+124Sn) reaction are larger than
those in the 40Ca+132Sn (48Ca+132Sn) reaction. The reason of that is the
nonzero quadrupole deformation of the heavy nucleus 124Sn. It should be
stressed that there are almost no difference between $\sigma_{cap}^{red}$ in
the reactions 40,48Ca+124,132Sn at energies above the Coulomb barrier.

Figure 2: The calculated capture cross sections versus $E_{\rm c.m.}$ for the
reactions 40Ca+124Sn (solid line) and 48Ca+124Sn (dashed line). The
experimental data for the reactions 40Ca+124Sn (solid squares) and 48Ca+124Sn
(open squares) are from Ref. Kol . In the calculations the barriers were
adjusted to the experimental values.  Figure 3: The calculated capture cross
sections versus $E_{\rm c.m.}$ for the reactions 40Ca+132Sn (solid line) and
48Ca+132Sn (dashed line). The experimental data for the reactions 40Ca+132Sn
(solid squares) and 48Ca+132Sn (open squares) are from Ref. Kol . In the
calculations the barriers were adjusted to the experimental values.

In Figs. 2 and 3 one can see a good agreement between the calculated results
and the experimental data in the reactions 40,48Ca+124,132Sn. This means that
the observed capture enhancements in the reactions 40Ca+124,132Sn at sub-
barrier energies are related to the two-neutron transfer effect. Note that the
slope of the excitation function strongly depends on the deformations of the
interacting nuclei and, respectively, on the neutron transfer effect.

To describe the reactions 40,48Ca+132Sn (Fig. 2) and 48Ca+124,132Sn (Fig. 3),
we extracted the values of the corresponding Coulomb barrier $V_{b}$ for the
spherical nuclei. There are differences between the calculated and extracted
$V_{b}$. From the direct calculations of the nucleus-nucleus potentials (with
the same set of parameters), we obtained
$V_{b}$(40Ca+124Sn)-$V_{b}$(48Ca+124Sn)=2.3 MeV,
$V_{b}$(40Ca+132Sn)-$V_{b}$(48Ca+132Sn)=2.2 MeV,
$V_{b}$(40Ca+124Sn)-$V_{b}$(40Ca+132Sn)=1.3 MeV, and
$V_{b}$(48Ca+124Sn)-$V_{b}$(48Ca+132Sn)=1.2 MeV. From the extractions, we got
$V_{b}$(40Ca+124Sn)-$V_{b}$(48Ca+124Sn)=1.1 MeV
$V_{b}$(40Ca+132Sn)-$V_{b}$(48Ca+132Sn)=1.0 MeV,
$V_{b}$(40Ca+124Sn)-$V_{b}$(40Ca+132Sn)=-0.3 MeV, and
$V_{b}$(48Ca+124Sn)-$V_{b}$(48Ca+132Sn)=-0.4 MeV, which seem to be
unrealistically small. However, these differences of $V_{b}$ do not influence
the slopes of the excitation functions but only lead to the shifting of the
energy scale. With realistic isospin trend of $V_{b}$
$\sigma_{cap}$(40Ca+124Sn)$<\sigma_{cap}$(48Ca+124Sn) and
$\sigma_{cap}$(40Ca+132Sn)$<\sigma_{cap}$(48Ca+132Sn) at energies above the
corresponding Coulomb barriers.

###  III.2 Effect of neutron transfer in reactions with beams 16,18O

Figures 4-7 show the capture excitation function for the reactions
16,18O+76,74Ge, 16,18O+94,92Mo, 16,18O+114,112,120,118,126,124Sn, and
16,18O+52,50Cr as a function of bombarding energy. One can see a rather good
agreement between the calculated results and the experimental data Jia ;
16OAGe ; AO92Mo ; AOASn  for the reactions 16O+76Ge, 16,18O+92Mo, and
18O+112,118,124Sn.

Figure 4: (Colour online) The calculated (solid line) capture cross sections
versus $E_{\rm c.m.}$ for the reactions 16O+76Ge and 18O+74Ge (the curves
coincide). For the 18O+74Ge reaction, the calculated capture cross sections
without neutron transfer are shown by dotted line. The experimental data for
the reactions 16O+76Ge (open circles) and 18O+74Ge (open squares) are from
Ref. Jia . The experimental data for the 16O+76Ge reaction (solid circles) are
from Ref. 16OAGe .  Figure 5: (Colour online) The calculated capture cross
sections versus $E_{\rm c.m.}$ for the reactions 16O+92Mo (dashed line) and
18O+92Mo (solid line). For the 18O+92Mo reaction, the calculated capture cross
sections without the neutron transfer are shown by dotted line. The
experimental data for the reactions 16O+92Mo (solid stars) and 18O+92Mo (solid
squares) are from Ref. AO92Mo .  Figure 6: The calculated capture cross
sections versus $E_{\rm c.m.}$ for the reactions 16O+114Sn and 18O+112Sn
(solid line), 16O+120Sn and 18O+118Sn (dashed line), 16O+126Sn and 18O+124Sn
(dotted line). The calculated results for the reactions 16O+114,120,126Sn and
18O+112,118,124Sn coincide, respectively. The experimental data for the
reactions 18O+112Sn (solid squares), 18O+118Sn (open squares), and 18O+124Sn
(open stars) are from Ref. AOASn .  Figure 7: (Coulor online) The calculated
capture cross sections versus $E_{\rm c.m.}$ for the reactions 16O+52Cr
(dashed line) and 18O+50Cr (solid line).

The $Q_{2n}$-values for the $2n$-transfer processes are positive (negative)
for all reactions with 18O (16O). Thus, the neutron transfer can be important
for the reactions with the 16O beam. However, our results show that cross
sections for reactions 16O+76Ge (16O+114,120,126Sn,52Cr) and 18O+74Ge
(18O+114,118,124Sn,50Cr) are very similar. The reason of such behavior is that
after the $2n$-transfer in the system 18O+A-2X$\to^{16}$O+AX the deformations
remain to be similar. As a result, the corresponding Coulomb barriers of the
systems 18O+A-2X and 16O+AX are almost the same and, correspondingly, their
capture cross sections coincide. Just the same behavior was observed in the
recent experiments 16,18O+76,74Ge Jia .

One can see in Figs. 4-7 that at energies above and near the Coulomb barrier
the cross sections with and without two-neutron transfer are quite similar.
After the $2n$-transfer (before the capture) in the reactions
18O($\beta_{2}=0.1$) + 92Mo($\beta_{2}=0.05$)$\to^{16}$O($\beta_{2}=0$) +
94Mo($\beta_{2}=0.151$), 18O($\beta_{2}=0.1$) +
74Ge($\beta_{2}=0.283$)$\to^{16}$O($\beta_{2}=0$) + 76Ge($\beta_{2}=0.262$),
18O($\beta_{2}=0.1$)+112Sn($\beta_{2}=0.123$)$\to^{16}$O($\beta_{2}=0$)+114Sn($\beta_{2}=0.121$),
18O($\beta_{2}=0.1$)+118Sn($\beta_{2}=0.111$)$\to^{16}$O($\beta_{2}=0$)+120Sn($\beta_{2}=0.104$),
and
18O($\beta_{2}=0.1$)+124Sn($\beta_{2}=0.095$)$\to^{16}$O($\beta_{2}=0$)+126Sn($\beta_{2}=0.09$)
the deformations of the nuclei decrease and the values of the corresponding
Coulomb barriers increase. As a result, the transfer suppresses the capture
process at the sub-barrier energies. The suppression becomes stronger with
decreasing energy. As examples, in Fig. 4 and 5 we show this effect for the
reactions 18O+74Ge,92Mo.

###  III.3 Capture and quasielastic barrier distributions

In Figs. 8 and 9, the calculated capture barrier distributions

$D=d^{2}(E_{\rm c.m.}\sigma_{cap})/dE_{\rm c.m.}^{2}$

for the reactions 16O+76Ge,144,154Sm have only one pronounced maximum around
$E_{\rm c.m.}=V_{b}$ as in the experiments Jia ; Timmers . The calculated
barrier distributions in Figs. 8 and 9 are slightly wider and fit the
experimental data better than those obtained with the couple-channels approach
in Fig. 5 of Ref. Jia . The capture (fusion) cross sections for the reactions
16O+76Ge,144,154Sm were well described with the quantum diffusion model in
Ref. EPJSub1 . With almost spherical (deformed) target-nucleus we obtain a
more narrow (wide) barrier distribution for the 16O+144Sm (16O+154Sm)
reaction.

We compared the capture and the quasielastic barrier distributions for these
reactions (Figs. 8 and 9).

Figure 8: (Colour online) (a) The calculated values of the quasielastic $\pi
R_{b}^{2}D_{qe}$ (solid line) and capture $D$ (dotted line) barrier
distributions for the reactions 16O + 76Ge and 18O + 74Ge. The curves coincide
for these reactions. The calculated $D$ for the spherical interacting nuclei
is shown by dashed line. The experimental data for the reactions 16O + 76Ge
(solid circles) and 18O + 74Ge (open circles) are from Ref. Jia . (b) The
calculated values of $\pi R_{b}^{2}D_{qe}$ (solid line) and $D$ (dotted line)
are shown in the logarithmic scale.  Figure 9: (Colour online) The calculated
values of quasielastic $D_{qe}$ (solid line) and capture $D/(\pi R_{b}^{2})$
(dotted line) barrier distributions for the reactions 16O + 144Sm (a) and 16O
+ 154Sm (b). The experimental $D_{qe}$ (open squares) and $D/(\pi R_{b}^{2})$
(solid circles) are from Ref. Timmers . The calculated $D$ for the spherical
interacting nuclei is shown by dashed line (b).

There is a direct relationship between the capture and the quasielastic
scattering processes because any loss from the quasielastic channel
contributes directly to the capture (the conservation of the reaction flux):

$P_{qe}(E_{\rm c.m.},J)+P_{cap}(E_{\rm c.m.},J)=1$

and

$dP_{cap}/dE_{\rm c.m.}=-dP_{qe}/dE_{\rm c.m.},$

where $P_{qe}$ is the reflection probability and $P_{cap}$ is the capture
(transmission) probability ($J$ is the partial wave). The quasielastic barrier
distribution is extracted by taking the first derivative of the $P_{qe}(E_{\rm
c.m.},J=0)$ or $P_{cap}(E_{\rm c.m.},J=0)$ with respect to $E_{\rm c.m.}$,
that is,

$D_{qe}(E_{\rm c.m.})=-dP_{qe}(E_{\rm c.m.},J=0)/dE_{\rm c.m.}=dP_{cap}(E_{\rm
c.m.},J=0)/dE_{\rm c.m.}.$

So, by employing the quantum diffusion approach and calculating
$dP_{cap}(E_{\rm c.m.},J=0)/dE_{\rm c.m.}$, one can obtain $D_{qe}(E_{\rm
c.m.})$. One can see in Figs. 8 and 9 that the shapes of the quasielastic and
capture barrier distributions are similar. The same conclusion was
experimentally obtained for the 20Ne+208Pb reaction in Ref. Piasecki . As in
the case of capture barrier distribution, one can show that the width of the
quasielastic barrier distribution increases with the deformation of the
target-nucleus. In addition to the mean peak position of the $D_{qe}$ around
the barrier height, we observe the sharp change of the slope of $D_{qe}$ or
$D$ below the Coulomb barrier energy because of a change of the regime of
interaction (the external turning point leaves the region of the nuclear
forces and friction EPJSub ; EPJSub1 ) in the deep sub-barrier capture process
(Fig. 8(b)).

##  IV Summary

As shown with the quantum diffusion approach, the capture cross sections for
the reactions 16O+52Cr,76Ge,94Mo,114,120,126Sn and
18O+50Cr,74Ge,92Mo,112,118,124Sn, respectively, almost match. The fusion
enhancement due to the positive $Q_{2n}$-value $2n$-transfer for 18O+74Ge has
not been observed Jia  because the deformations of nuclei slightly decrease
after the neutron transfer. This is different from the situation for the
reactions 40Ca+124,132Sn Kol  with large positive $Q_{2n}$-values. The strong
enhancements have been observed Kol  in these reactions at sub-barrier
energies because the deformation of light nucleus strongly increases (the
heavy nucleus is spherical before and after transfer) after the two-neutron
transfer.

We found that the shapes of the quasielastic and capture barrier distributions
are similar. The sharp change of the slope of the quasielastic or capture
barrier distribution is predicted at deep sub-barrier energy. This anomalous
behavior of the barrier distribution is expected to be the experimental
indication of a change of the regime of interaction in the sub-barrier
capture. One concludes that the quasielastic technique could be an important
tool in capture (fusion) research.

This work was supported by DFG, NSFC, and RFBR. The IN2P3(France) -
JINR(Dubna) and Polish - JINR(Dubna) Cooperation Programmes are gratefully
acknowledged. H.Q. Zhang is grateful to Chinese NSFC for the partial support.

## References

  * (1) G. Montagnoli, et al., Phys. Rev. C 85, 024607 (2012); C. Simenel, Eur. Phys. J. A 48, 152 (2012); C.A. Bertulani, EPJ Web Conf. 17, 15001 (2011); Z. Kohley et al., Phys. Rev. Lett. 107, 202701 (2011); J.F. Liang, EPJ Web Conf. 17, 02002 (2011); F. Scarlassara et al., EPJ Web Conf. 17, 05002 (2011).
  * (2) H.M. Jia et al., Phys. Rev. C 86, 044621 (2012).
  * (3) J.J. Kolata et al., Phys. Rev. C 85, 054603 (2012).
  * (4) V.V. Sargsyan, G.G. Adamian, N.V. Antonenko, and W. Scheid, Eur. Phys. J. A 45, 125 (2010); Eur. Phys. J. A 47, 38 (2011); Eur. Phys. J. A 48, 118 (2012); Eur. Phys. J. A 49, 19 (2013).
  * (5) V.V. Sargsyan, G.G. Adamian, N.V. Antonenko, W. Scheid, and H.Q. Zhang, Phys. Rev. C 84, 064614 (2011); Phys. Rev. C 85, 024616 (2012); Phys. Rev. C 86, 014602 (2012).
  * (6) S. Ayik, B. Yilmaz, and D. Lacroix, Phys. Rev. C 81, 034605 (2010).
  * (7) H. Hofmann, Phys. Rep. 284, 137 (1997); S. Ayik, B. Yilmaz, A. Gokalp, O. Yilmaz, and N. Takigawa, Phys. Rev. C 71, 054611 (2005); V.V. Sargsyan , Z. Kanokov, G.G. Adamian, and N.V. Antonenko, Part. Nucl. 41, 175 (2010); G. Hupin and D. Lacroix, Phys. Rev. C 81, 014609 (2010).
  * (8) R.A. Broglia, C.H. Dasso, S. Landowne, and A. Winther, Phys. Rev. C 27, 2433 (1983); R.A. Broglia, C.H. Dasso, S. Landowne, and G. Pollarolo, Phys. Lett. B 133, 34 (1983).
  * (9) S. Szilner et al., Phys. Rev. C 76, 024604 (2007); S. Szilner et al., Phys. Rev. C 84, 014325 (2011); L. Corradi et al., Phys. Rev. C 84, 034603 (2011).
  * (10) G.G. Adamian, A.K. Nasirov, N.V. Antonenko, and R.V. Jolos, Phys. Part. Nucl. 25, 583 (1994); K. Washiyama, D. Lacroix, and S. Ayik, Phys. Rev. C 79, 024609 (2009); S. Ayik, K. Washiyama, and D. Lacroix, Phys. Rev. C 79, 054606 (2009).
  * (11) S. Raman, C.W. Nestor, Jr, and P. Tikkanen, At. Data Nucl. Data Tables 78, 1 (2001).
  * (12) L.F. Canto, P.R.S. Gomes, J. Lubian, L.C. Chamon, and E. Crema, J. Phys. G 36, 015109 (2009); Nucl. Phys. A821, 51 (2009).
  * (13) E.F. Aguilera, J.J. Kolata, and R.J. Tighe, Phys. Rev. C 52, 3103 (1995).
  * (14) M. Benjelloun, W. Galster, and J. Vervier, Nucl. Phys. A560, 715 (1993).
  * (15) P. Jacobs, Z. Fraenkel, G. Mamane, and L. Tserruya, Phys. Lett. B 175, 271 (1986).
  * (16) H. Timmers et al., Nucl. Phys. A584, 190 (1995).
  * (17) E. Piasecki et al., Phys. Rev. C 85, 054608 (2012).