# Input Delay Compensation for Neuron Growth by PDE Backstepping

Cenk Demir  Shumon Koga  Miroslav Krstic  University of California, San Diego,
CA 92093 USA, (e-mail: cdemir@ucsd.edu).  University of California, San Diego,
CA 92093 USA, (e-mail: skoga@ucsd.edu).  University of California, San Diego,
CA 92093 USA, (e-mail: krstic@ucsd.edu).

###### Abstract

Neurological studies show that injured neurons can regain their functionality
with therapeutics such as Chondroitinase ABC (ChABC). These therapeutics
promote axon elongation by manipulating the injured neuron and its
intercellular space to modify tubulin protein concentration. This fundamental
protein is the source of axon elongation and its spatial distribution is the
state of the axon growth dynamics. Such dynamics often contain time delays
because of biological processes. This work introduces an input delay
compensation with state-feedback control law for axon elongation by regulating
tubulin concentration. Axon growth dynamics with input delay is modelled as
coupled parabolic diffusion-reaction-advection Partial Differential Equations
(PDE) with a boundary governed by Ordinary Differential Equations (ODE),
associated with a transport PDE. A novel feedback law is proposed by using
backstepping method for input-delay compensation. The gain kernels are
provided after transforming the interconnected PDE-ODE-PDE system to a target
system. The stability analysis is presented by applying Lyapunov analysis to
the target system in the spatial $\mathcal{H}_{1}$-norm, thereby the local
exponential stability of the original error system is proved by using norm
equivalence.

###### keywords:

Backstepping, input delay, moving boundary, axon elongation, Stefan problem

##  1 Introduction

The study about axonal growth is an expanding field in neuroscience in order
to understand functionality and structure of neurons (Yamada et al., 1970;
Tessier-Lavigne and Goodman, 1996; Kandel et al., 2000). This understanding
helps to cure neurological disorders such as spinal cord injury and
Alzheimer’s disease (Liu et al., 1997; Maccioni et al., 2001). For example,
Chondroitinase ABC is a therapeutic in clinical trials to cure spinal cord
injuries. It aims to restore axon growth in damaged neurons by manipulating
the extracellular matrix, the non-cellular macromolecules surrounding the cell
(Karimi-Abdolrezaee et al., 2010; Bradbury and Carter, 2011; Frantz et al.,
2010).

Neurons are specialized cells in the nervous system that receive and transmit
electrical signals. These signals generate perception, and initiate physical
actions. Neurons consist of three main components that facilitate perception
and give commands to the muscles; the soma, the axon and the growth cone as
shown in Fig. 1. The soma is the body of the neuron and is responsible for
producing proteins. Dendrites, which resemble tree branches, extend from the
soma and receive incoming electrical signals (Kandel et al., 2000). The axon
is a wire-like structure that electrical signals travel along. The growth cone
is the highly mobile structure at the tip of the axon that seeks chemical cues
from the target neuron which will receive the signal. The signal transmission
process starts with the entry of an electric signal from the dendrites. The
signal then moves along the axon to the growth cone. In the final step, the
growth cone identifies the postsynaptic neuron that will receive the signal
(Julien, 1999). With the transmission of the signal, the process is complete.

Tubulin proteins play an essential role in signal transmission by serving as
the building blocks which extend the axon. Tubulin monomers and dimers build
chemical bonds to create microtubules, the tubular protein shape which forms
the structure of the axon (Desai and Mitchison, 1997). The creation of
microtubules depends on the dynamics of tubulin concentration in the neuron.
These dynamics include tubulin production rate, assembly-disassembly rate, and
the tubulin transportation process (Diehl et al., 2014). Neurological
disorders and external damage, such as spinal cord injuries, can increase the
number of axon growth inhibitors that prevent axon elongation (Lemons et al.,
1999). In ChABC therapy, bacterial enzymes are injected into damaged tissue to
digest axon growth inhibitors which prevent the formation of microtubules (Lee
et al., 2010; Frantz et al., 2010). The potential to control microtubule
formation motivated researchers to design a control law that manipulates axon
elongation to achieve the desired lengths (Demir et al., 2021). Although Demir
et al. (2021) provides feedback control law, it does not account for the input
time delay between the injected input and the tubulin concentration dynamics.
Time delays frequently occur in biological processes and increase the
complexity of the system. In this paper, we deal with this problem by applying
an input delay compensation with a novel feedback control law.

There are several mathematical models which seek to describe the behavior of
tubulin concentration. Van Veen and Van Pelt (1994) presents a nonlinear ODE
model which includes the tubulin production at the soma, tubulin
transportation, microtubule formation, and axon elongation. Another model
expresses the microtubule formation process by linking the mechanics of the
membrane and the axon growth (García-Grajales et al., 2017). Besides these ODE
models, researchers also developed PDE models that enhance the understanding
of the behavior of tubulin concentration and axon elongation (McLean et al.,
2004; Diehl et al., 2014). In this model, the authors use a PDE model to
describe tubulin concentration along the axon and a nonlinear ODE to describe
the axon length. In addition, there is a coupling between PDE and ODE because
the length of the axon is the domain of the tubulin concentration in this
model. Moreover, McLean and Graham (2006) provides a stability analysis of
this coupled PDE-ODE system.

Over the last two decades, the boundary control strategies of PDE-based models
have been extensively studied for parabolic, hyperbolic, and some other exotic
PDEs (Krstic and Smyshlyaev, 2008). The contribution of Smyshlyaev and Krstic
(2004) allowed analytically unsolvable kernel PDEs to be solved numerically by
the method of successive approximation, which commenced the usage of boundary
control for relatively complex problems. Then, (Susto and Krstic, 2010; Tang
and Xie, 2011) extended the boundary control law for the class of coupled PDE
and ODE systems. Following this extension, (Krstic, 2009b) introduced input
delay compensation for boundary control by considering an input delay as a
transport PDE. Krstic (2009a) also provides input delay compensation and
control for unstable reaction-diffusion PDE for arbitrarily long input delays.
Building on previous studies which examined PDEs with a constant domain size
in time, Koga et al. (2018) and Koga and Krstic (2020) developed a
backstepping-based boundary control law with global stability results for
parabolic PDEs with the moving boundary, called Stefan problem. In addition,
delay compensation for the one phase Stefan problem presented in Koga et al.
(2020). Buisson-Fenet et al. (2018); Yu et al. (2020) obtained local stability
for the class of coupled hyperbolic PDE with a moving boundary governed by an
ODE for a piston movement and traffic congestion problems, respectively. In
addition, Demir et al. (2021, 2022) provide the local stability results for a
parabolic PDE with a moving boundary governed by a nonlinear ODE for the axon
growth problem.

This paper presents the feedback stabilization with input delay compensation
for the coupled model consisting of tubulin concentration and axon growth.
First, we introduce the input delay as a transport PDE. We obtain the
reference error system by subtracting the steady-state solution from the
dynamics and use linearization to ODE state around the zero state to deal with
nonlinearity. Then, we apply backstepping transformations to the transport PDE
and the parabolic PDE. While we derive some of the kernels of these
transformations analytically, we obtain some of them numerically. By using
these kernel solutions, we determine the control law. Finally, we prove the
local stability of the original PDE-ODE-PDE system.

This paper is structured as follows. Section II presents the modeling of
tubulin concentration and axon growth as a coupled PDE-ODE model with input
delay. Section III presents the feedback control law with compensating input
delay by using the method of backstepping. The next section proves the
stability of the closed-loop system. The paper ends with the conclusion in
Section IV.

Figure 1: Schematic of neuron and state variables

##  2 Modeling and Problem Statement

###  2.1 Axon Growth Model

In this model, tubulin protein is considered as a responsible factor for axon
growth with two assumptions. The assumptions are that tubulin is modeled as a
homogeneous continuum and that only tubulin is effective for axon growth.
Then, axon growth and tubulin concentration with input delay is modeled as

$\displaystyle c_{t}(x,t)=$ $\displaystyle Dc_{xx}(x,t)-ac_{x}(x,t)-gc(x,t),$
(1) $\displaystyle c_{x}(0,t)=$ $\displaystyle-q_{\rm s}(t-D_{e}),$ (2)
$\displaystyle c(l(t),t)=$ $\displaystyle c_{\rm c}(t),$ (3) $\displaystyle
l_{\rm c}\dot{c}_{\rm c}(t)=$ $\displaystyle(a-gl_{\rm c})c_{\rm
c}(t)-Dc_{x}(l(t),t)$ $\displaystyle-(r_{\rm g}c_{\rm c}(t)+\tilde{r}_{\rm
g}l_{\rm c})(c_{\rm c}(t)-c_{\infty}),$ (4) $\displaystyle\dot{l}(t)=$
$\displaystyle r_{\rm g}(c_{c}(t)-c_{\infty}),$ (5)

In (1)-(5), $c(x,t)$ represents tubulin concentration along the axon, $q_{s}$
represents the tubulin flux and subscripts $s$, and $c$ are denoted for soma
and growth cone respectively. Namely, $c_{c}(t)$ denotes the tubulin flux in
the growth cone, and $q_{s}(t)$ denotes tubulin flux in soma. The axon length
is described as $l(t)$ which is in $x$-coordinate. The constants $D$, $a$, and
$g$, denote the tubulin diffusion constant, velocity constant, and degradation
rate, respectively. $l_{\rm c}$ represents the growth ratio, $r_{\rm g}$ is a
lumped parameter, and $\tilde{r}_{\rm g}$ is the microtubule reaction rate.
$c_{\infty}$ is the equilibrium of the tubulin concentration in the growth
cone. The last constant $D_{e}$ represents the time delay of input.

###  2.2 Problem Statement

We pursue to drive the axon length to a given desired length $l_{s}>0$ by
designing the state feedback control law of $q_{\rm s}(t-D_{e})$ for known
time delay, $D_{e}$. The steady-state solution of the system (1)–(5) is then
derived by considering an equilibrium of the tubulin concentration $c_{\rm
eq}(x)$. Thus, the control objective is to achieve

$\displaystyle\lim_{t\to\infty}l(t)$ $\displaystyle=l_{\rm s},$ (6)
$\displaystyle\lim_{t\to\infty}c(x,t)$ $\displaystyle=c_{eq}(x),$ (7)

through compensating the input delay $D_{e}$.

###  2.3 Steady-State Solution

The steady-state solution of the system (1)-(5) is solved for the desired axon
length $l_{\rm s}$, which is, as derived in Demir et al. (2021), given by

$\displaystyle c_{\rm eq}(x)=c_{\infty}\left(K_{+}e^{\lambda_{+}(x-l_{\rm
s})}+K_{-}e^{\lambda_{-}(x-l_{\rm s})}\right),$ (8) $\displaystyle q_{\rm
s}^{*}=-c_{\infty}\left(K_{+}\lambda_{+}e^{-\lambda_{+}l_{\rm
s}}+K_{-}\lambda_{-}e^{-\lambda_{-}l_{\rm s}}\right),$ (9)

where

$\displaystyle\lambda_{+}=$
$\displaystyle\frac{a+\sqrt{a^{2}+4Dg}}{2D},\quad\lambda_{-}=\frac{a-\sqrt{a^{2}+4Dg}}{2D},$
(10) $\displaystyle K_{+}=$ $\displaystyle\frac{1}{2}+\frac{a-2gl_{\rm
c}}{2\sqrt{a^{2}+4Dg}},\ K_{-}=\frac{1}{2}-\frac{a-2gl_{\rm
c}}{2\sqrt{a^{2}+4Dg}}.$ (11)

###  2.4 Reference Error System

The reference error states, denoted as $u(x,t)$, $z_{1}(t)$, $z_{2}(t)$ and
$U(t-D_{e})$, are defined by

$\displaystyle u(x,t)=$ $\displaystyle c(x,t)-c_{\rm eq}(x),$ (12)
$\displaystyle z_{1}(t)=$ $\displaystyle c_{\rm c}(t)-c_{\infty},$ (13)
$\displaystyle z_{2}(t)=$ $\displaystyle l(t)-l_{\rm s},$ (14) $\displaystyle
U(t-D_{e})=$ $\displaystyle-(q_{\rm s}(t-D_{e})-q_{\rm s}^{*}).$ (15)

From (1)-(5), by using (8)-(9), and (12)-(15), we obtain

$\displaystyle u_{t}(x,t)=$ $\displaystyle Du_{xx}(x,t)-au_{x}(x,t)-gu(x,t),$
(16) $\displaystyle u_{x}(0,t)=$ $\displaystyle U(t-D_{e}),$ (17)
$\displaystyle u(l(t),t)=$ $\displaystyle z_{1}(t)+\tilde{h}(z_{2}(t)),$ (18)
$\displaystyle\dot{z}_{1}(t)=$ $\displaystyle\tilde{a}z_{1}(t)-\beta
u_{x}(l(t),t)-\kappa z_{1}(t)^{2},$ (19) $\displaystyle\dot{z}_{2}(t)=$
$\displaystyle r_{\rm g}z_{1}(t),$ (20)

where

$\displaystyle\tilde{a}=\frac{a-r_{\rm g}c_{\infty}}{l_{\rm
c}}-g-\tilde{r}_{\rm g},\quad\beta=\frac{D}{l_{\rm
c}},\quad\kappa=\frac{r_{\rm g}}{l_{\rm c}},$ (21)
$\displaystyle\tilde{h}(z_{2}(t))=c_{\infty}\left(1-K_{+}e^{\lambda_{+}z_{2}(t)}-K_{-}e^{\lambda_{-}z_{2}(t)}\right).$
(22)

The delayed controller can be represented by the following transport PDE
dynamics:

$\displaystyle v_{t}(x,t)=$ $\displaystyle v_{x}(x,t),\quad D_{e}\geq x\geq 0$
(23) $\displaystyle v(D_{e},t)=$ $\displaystyle U(t).$ (24)

The solution to this equation is $v(x,t)=U(t+x-D_{e})$ and thus the output
$v(0,t)=U(t-D_{e})$ gives the delayed input. Then, let $X\in\mathbb{R}^{2}$ be
an ODE state vector for the reference error states $z_{1}(t)$ and $z_{2}(t)$,
defined by

$\displaystyle X(t)=[z_{1}(t)\quad z_{2}(t)]^{\top}.$ (25)

By using (25), the systematic representation of (16)-(20) and (23)-(24) can be
written as

$\displaystyle u_{t}(x,t)=$ $\displaystyle Du_{xx}(x,t)-au_{x}(x,t)-gu(x,t),$
(26) $\displaystyle u_{x}(0,t)=$ $\displaystyle v(0,t),$ (27) $\displaystyle
u(l(t),t)=$ $\displaystyle e_{1}X(t)+\tilde{h}(e_{2}X(t)),$ (28)
$\displaystyle\dot{X}(t)=$ $\displaystyle AX(t)+f(X(t))+Bu_{x}(l(t),t),$ (29)
$\displaystyle v_{t}(x,t)=$ $\displaystyle v_{x}(x,t),\qquad\qquad D_{e}\geq
x\geq 0$ (30) $\displaystyle v(D_{e},t)=$ $\displaystyle U(t).$ (31)

where

$\displaystyle A=\left[\begin{array}[]{cc}\tilde{a}&0\\\ r_{\rm
g}&0\end{array}\right],\ B=\left[\begin{array}[]{c}-\beta\\\
0\end{array}\right],$ (36) $\displaystyle
f(X(t))=-\kappa(e_{1}X(t))^{2}e_{1},$ (37)

where $e_{1}$ and $e_{2}$ are unit vectors; $e_{1}=[1\quad 0]$, $e_{2}=[0\quad
1]$.

###  2.5 Linearized Reference Error System

Applying linearization around zero states to (26)-(31) leads us to the
following linearized system:

$\displaystyle\dot{X}(t)=$ $\displaystyle AX(t)+Bu_{x}(l(t),t),$ (38)
$\displaystyle u_{t}(x,t)=$ $\displaystyle Du_{xx}(x,t)-au_{x}(x,t)-gu(x,t),$
(39) $\displaystyle u(l(t),t)=$ $\displaystyle C^{\top}X(t),$ (40)
$\displaystyle u_{x}(0,t)=$ $\displaystyle v(0,t),$ (41) $\displaystyle
v_{t}(x,t)=$ $\displaystyle v_{x}(x,t),\quad D_{e}\geq x\geq 0$ (42)
$\displaystyle v(D_{e},t)=$ $\displaystyle U(t).$ (43)

where

$\displaystyle C=\left[1\quad-\frac{(a-gl_{\rm
c})c_{\infty}}{D}\right]^{\top}.$ (44)

##  3 State Feedback Design

Before we start introducing the state feedback control law, we use the
following notations for several norms:

$\displaystyle||u(\cdot,t)||_{L_{2}}=\sqrt{\int_{0}^{l(t)}u(.,t)^{2}dx},$
$\displaystyle||u(\cdot,t)||_{H_{1}}=\sqrt{\int_{0}^{l(t)}u(.,t)^{2}+u_{x}(.,t)^{2}dx}.$

In the rest of the paper, we propose the backstepping control law that deals
with input delay compensation and stability proof.

###  3.1 Backstepping Transformation

We consider backstepping transformation of the form

$\displaystyle w(x,t)=$ $\displaystyle u(x,t)-\int_{x}^{l(t)}k(x,y)u(y,t)dy$
$\displaystyle-\phi(x-l(t))X(t),$ (45) $\displaystyle z(x,t)=$ $\displaystyle
v(x,t)+\int_{0}^{x}p(x,y)v(y,t)dy$
$\displaystyle+\int_{0}^{l(t)}q(x,y)u(y,t)dy+\psi(x-l(t))X(t),$ (46)

where the kernels $k(x,y)$, $\phi(x)$, $p(x,y)$, $q(x,y)$ and $\psi(x)$ are to
be determined to transform the original system to a target system.

$\displaystyle\dot{X}(t)=$ $\displaystyle(A+BK)X(t)+Bw_{x}(l(t),t)+f(X(t)),$
(47) $\displaystyle w_{t}(x,t)=$ $\displaystyle
Dw_{xx}(x,t)-aw_{x}(x,t)-gw(x,t)$
$\displaystyle-\dot{l}(t)F(x,X(t))-\phi(x-l(t))^{\top}f(X(t))$
$\displaystyle-\left(\phi^{\prime}(x-l(t))^{\top}+\frac{a}{D}\phi(x-l(t))^{\top}\right)B\bar{h}(X),$
(48) $\displaystyle w(l(t),t)=$ $\displaystyle\bar{h}(X(t)),$ (49)
$\displaystyle w_{x}(0,t)=$ $\displaystyle z(0,t),$ (50) $\displaystyle
z_{t}(x,t)=$ $\displaystyle z_{x}(x,t),$
$\displaystyle-\dot{l}(t)\left(q(x,l(t))C^{\top}+\psi^{\prime}(x-l(t))\right)X(t),$
(51) $\displaystyle z(D_{e},t)=$ $\displaystyle 0,$ (52)

where the nonlinear terms are

$\displaystyle\tilde{h}(z_{2}(t))=$ $\displaystyle
c_{\infty}\left(1-K_{+}e^{\lambda_{+}z_{2}(t)}-K_{-}e^{\lambda_{-}z_{2}(t)}\right),$
(53) $\displaystyle\bar{h}(X)=$
$\displaystyle\left(z_{1}(t)+\tilde{h}(z_{2}(t))\right)-C^{\top}X,$ (54)

where $K\in\mathbb{R}^{2}$ is a chosen feedback control gain vector to make
$A+BK$ Hurwitz. Thus, $K$ is

$\displaystyle K=[k_{1}\quad k_{2}]^{\top},\quad
k_{1}>\frac{\tilde{a}}{\beta},\quad k_{2}>0.$ (55)

The redundant nonlinear term is

$\displaystyle
F(x,X(t))=\left(\phi^{\prime}(x-l(t))^{T}-k(x,l(t))C^{T}\right)X(t).$ (56)

###  3.2 Gain Kernel Solution

The conditions of the kernel functions are obtained to satisfy both the
governing equations. Thus, we take the time and spatial derivatives of (45)
associated to the solutions of (38)-(41), and then we have the following gain
kernel solutions [as detailed in Demir et al. (2021)]

$\displaystyle\phi(x)^{\top}=\begin{bmatrix}C^{\top}&K^{\top}-\frac{1}{D}C^{\top}BC^{\top}\end{bmatrix}e^{N_{1}x}\begin{bmatrix}I\\\
0\end{bmatrix},$ (57)

where the matrix $N_{1}\in\mathbb{R}^{4\times 4}$ is defined as

$\displaystyle
N_{1}=\begin{bmatrix}0&\frac{1}{D}\left(gI+A+\frac{a}{D}BC^{\top}\right)\\\
I&\frac{1}{D}\left(BC^{\top}+aI\right)\end{bmatrix},$ (58)

and

$\displaystyle k(x,y)=\frac{1}{D}\phi(x-y)^{\top}B.$ (59)

Next, we derive the kernel functions in another transformation (46). Taking
time and spatial derivatives of (46) along with (42)-(43), we have the
following coupled PDE-ODE for gain kernels

$\displaystyle q_{x}(x,y)=$ $\displaystyle Dq_{yy}(x,y)+aq_{y}(x,y)-gq(x,y),$
(60) $\displaystyle q(0,y)=$ $\displaystyle-k_{x}(0,y),$ (61) $\displaystyle
q(x,l(t))=$ $\displaystyle-\frac{1}{D}\psi(x-l(t))B,$ (62) $\displaystyle
q_{y}(x,0)=$ $\displaystyle-\frac{a}{D}q(x,0),$ (63)
$\displaystyle\psi^{\prime}(x-l(t))=$
$\displaystyle\psi(x-l(t))A-Dq_{y}(x,l(t))C^{\top}$ $\displaystyle-
aq(x,l(t))C^{\top},$ (64) $\displaystyle\psi(-l(t))=$
$\displaystyle-\phi^{\prime}(-l(t)).$ (65)

Through the variable change, and the method of successive approximation, we
can show that there exists a unique classic solution for these equations.
Then, the last kernel function, $p(x,y)$, is the solution to the following
trasport PDE

$\displaystyle p_{x}(x,y)=$ $\displaystyle-p_{y}(x,y),$ (66) $\displaystyle
p(x,0)=$ $\displaystyle-Dq(x,0).$ (67)

###  3.3 Inverse Transformation

The inverse transformation of (45)-(46) is formulated as

$\displaystyle u(x,t)$
$\displaystyle=w(x,t)+\int_{x}^{l(t)}\iota(x,y)w(y,t)dy$ $\displaystyle\ \
+\theta(x-l(t))X(t),$ (68) $\displaystyle v(x,t)$
$\displaystyle=z(x,t)-\int_{0}^{x}\varrho(x,y)z(y,t)dy$ $\displaystyle\ \
-\int_{0}^{l(t)}\chi(x,y)w(y,t)dy-\varphi(x-l(t))X(t)$ (69)

Then, by a similar procedure to the derivation in the previous section, we
have the kernel functions $\iota(x,y)$, $\theta(x)$, $\varrho(x,y)$,
$\chi(x,y)$ and $\varphi(x)$ which are determined as the solution to the
following set of PDE-ODE equations:

$\displaystyle\iota_{xx}(x,y)-\iota_{yy}(x,y)=\frac{a}{D}\left(\iota_{x}(x,y)+\iota_{y}(x,y)\right),$
(70) $\displaystyle\iota_{x}(x,x)+\iota_{y}(y,y)=0,$ (71)
$\displaystyle\iota(x,l(t))=-\frac{1}{D}\theta(x-l(t))^{\top}B,$ (72)
$\displaystyle
D\theta^{{}^{\prime\prime}}(x-l(t))^{\top}+a\theta^{{}^{\prime}}(x-l(t))^{\top}$
$\displaystyle\qquad\qquad+\left(gI+A+BK^{\top}\right)\theta(x-l(t))^{\top}=0,$
(73) $\displaystyle\theta(0)=C,\quad\theta^{\prime}(0)=K,$ (74)

and

$\displaystyle\chi_{x}(x,y)$
$\displaystyle=-D\chi_{yy}(x,y)-a\chi_{y}(x,y)-g\chi(x,y),$ (75)
$\displaystyle\chi_{y}(x,o)$ $\displaystyle=-\frac{a}{D}\chi(x,0),$ (76)
$\displaystyle\chi(x,l(t))$
$\displaystyle=\frac{1}{D}\varphi(x-l(t))^{\top}B,$ (77)
$\displaystyle\chi(0,y)$ $\displaystyle=-\iota_{x}(0,y),$ (78)
$\displaystyle\varphi^{\prime}(x-l(t))^{\top}$
$\displaystyle=\varphi(x-l(t))^{\top}\left(A+BK^{\top}\right),$ (79)
$\displaystyle\varphi(-l(t))^{\top}$
$\displaystyle=\theta^{\prime}(-l(t))^{\top},$ (80)
$\displaystyle\varrho_{x}(x,y)$ $\displaystyle=\varrho_{y}(x,y),$ (81)
$\displaystyle\varrho(0,t)$ $\displaystyle=-\chi(x,0).$ (82)

Notice that (70)-(74) are analytically solvable, and the solution of these PDE
and ODE equations are detailed in Demir et al. (2021). In addition, (75)-(82)
are also analytically solvable, and the solutions are bounded.

###  3.4 Control Law

Finally, we consider the closed loop system consisting of the plant (47)-(52)
and the control law

$\displaystyle U(t)=$
$\displaystyle-\int_{t-D_{e}}^{t}p(D_{e},\vartheta+D_{e}-t)U(\vartheta)d\vartheta$
$\displaystyle-\int_{0}^{l(t)}q(D_{e},y)u(y,t)dy-\psi(D_{e}-l(t))X(t).$ (83)

(83) is obtained by using the boundary condition (52) in (46), and defining
$\vartheta=t+y-D_{e}$ where $\vartheta\in(t-D_{e},t)$ and
$v(y,t)=U(t+y-D_{e})$.

##  4 Stability Analysis

In this section, a local stability of the closed-loop system is presented by
considering the $H_{1}$-norm

$\displaystyle
Z(t)=||u(\cdot,t)||_{H_{1}}^{2}+||v(\cdot,t)||_{H_{1}}^{2}+X^{\top}X.$ (84)

Before presenting our major theorem, we impose the following assumptions.

###### Assumption 1

The axon length $l(t)$ maintains positive and is upper bounded, i.e., there
exists a positive constant $\bar{l}>0$ such that the following inequality
holds:

$\displaystyle 0<l(t)\leq\bar{l},$ (85)

for all $t\geq 0$.

###### Assumption 2

The time derivative of the axon length is also bounded, i.e., there exists a
positive constant $\bar{v}>0$ such that the following inequality holds:

$\displaystyle\big{|}\dot{l}(t)\big{|}\leq\bar{v},$ (86)

for all $t\geq 0$.

###### Assumption 3

The physical constants satisfy the following inequality:

$\displaystyle\frac{D}{4\bar{l}}\geq a.$ (87)

Assumption 3 is physically valid and applicable when the axon length is short
which is the case for the most of the neurons because average axon length is
$20\mu$m-4$0\mu$m (Clark et al., 2009). Then, we state our main result in the
following theorem.

###### Theorem 1

Let Assumption 3 hold. Consider the closed-loop system consisting of the plant
(26)–(31) with the control law (83). Then, there exist positive parameters
$\bar{K}>0$, $c>0$, and $\kappa>0$, such that if $Z(0)<\bar{K}$ then the
following norm estimate holds

$\displaystyle Z(t)\leq cZ(0)\exp(-\kappa t),$ (88)

for all $t\geq 0$, which guarantees the local exponential stability of the
closed-loop system.

To prove Theorem 1, first we apply the following transformation to (47)-(52)

$\displaystyle\varpi(x,t)=w(x,t)-xz(0,t)+l(t)z(0,t)-\bar{h}(X),$ (89)

so we have

$\displaystyle\dot{X}(t)=$
$\displaystyle(A+BK)X(t)+B\varpi_{x}(l(t),t)+Bz(0,t)$ $\displaystyle+f(X(t)),$
(90) $\displaystyle\varpi_{t}(x,t)=$ $\displaystyle
D\varpi_{xx}(x,t)-a\varpi_{x}(x,t)-g\varpi(x,t)$
$\displaystyle-\bar{h}^{\prime}(X)-g\bar{h}(X)-az(0,t)-gxz(0,t)$
$\displaystyle+gl(t)z(0,t)-xz_{t}(0,t)+\dot{l}(t)z(0,t)$
$\displaystyle+l(t)z_{t}(0,t)-\dot{l}(t)F(x,X(t))$
$\displaystyle-\phi(x-l(t))^{\top}f(X(t))$
$\displaystyle-\left(\phi^{\prime}(x-l(t))^{\top}B+\frac{a}{D}\phi(x-l(t))^{\top}B\right)\bar{h}(X),$
(91) $\displaystyle\varpi(l(t),t)=$ $\displaystyle 0,$ (92)
$\displaystyle\varpi_{x}(0,t)=$ $\displaystyle 0,$ (93) $\displaystyle
z_{t}(x,t)=$ $\displaystyle z_{x}(x,t)$
$\displaystyle-\dot{l}(t)\left(\psi^{\prime}(x-l(t))+q(x,l(t))C^{\top}\right)X(t),$
(94) $\displaystyle z(D_{e},t)=$ $\displaystyle 0,$ (95)

where $x\in(0,l(t))$ for (91)-(93), and $x\in[0,D_{e})$ for (94)-(95).

###  4.1 Proof of Lyapunov stability

We consider the following Lyapunov function for the target system

$\displaystyle V=d_{1}V_{1}+V_{2}+d_{2}V_{3}+d_{3}V_{4}+d_{4}V_{5},$ (96)

where $d_{1}>0$, $d_{2}>0$, $d_{3}>0$ and $d_{4}>0$, and each Lyapunov
function is

$\displaystyle V_{1}:=\frac{1}{2}$ $\displaystyle||\varpi||_{L_{2}}^{2},\quad
V_{2}:=\frac{1}{2}||\varpi_{x}||_{L_{2}}^{2},$ (97) $\displaystyle V_{3}:=$
$\displaystyle\frac{1}{2}\int_{0}^{D_{e}}e^{cx}z^{2}(x,t)dx,$ (98)
$\displaystyle V_{4}:=$
$\displaystyle\frac{1}{2}\int_{0}^{D_{e}}e^{cx}z_{x}^{2}(x,t)dx,$ (99)
$\displaystyle V_{5}:=$ $\displaystyle X^{\top}PX,$ (100)

where $c>0$, and $P>0$ is a positive definite matrix satisfying the Lyapunov
equation:

$\displaystyle(A+BK)^{\top}P+P(A+BK)=-Q,$ (101)

for some positive definite matrix $Q>0$. Since $A+BK^{\top}$ is Hurwitz, due
to the positive definiteness of $P$ and $Q$,

$\displaystyle\lambda_{\rm min}(P)X^{\top}X\leq X^{\top}PX\leq\lambda_{\rm
max}(P)X^{\top}X,$ (102)

holds, where $\lambda_{\max}(P)>\lambda_{\min}(P)>0$ are the largest and
smallest eigenvalues in the positive definite matrix $P$. Then, we state the
following lemma.

###### Lemma 1

Let Assumptions 1-2 and 3 are satisfied with

$\displaystyle\bar{v}=\frac{D-2\epsilon}{8\bar{l}},$ (103)

where $\epsilon$ can be picked $D>2\epsilon$. Then, for sufficiently large
$d_{1}>0$, $d_{2}>0$, $d_{3}>0$ and small $d_{4}>0$, there exist positive
constants $\beta_{1}>0$, and $\beta_{2}>0$ such the following norm estimate
holds for all $t\geq 0$:

$\displaystyle\dot{V}\leq-\alpha V+\beta_{1}V^{3/2}+\beta_{2}V^{2},$ (104)

where

$\displaystyle\alpha=$
$\displaystyle\min\left\\{d_{1}\frac{D}{8},\frac{g}{16},\frac{\lambda_{\rm
min}(Q)}{2\lambda_{\min}(P)},c\right\\},$ (105) $\displaystyle\beta_{1}=$
$\displaystyle d_{1}\frac{\bar{l}r_{\rm g}}{2}+d_{1}\frac{1}{2}r_{\rm
g}+d_{2}e^{cD_{e}}r_{\rm g}+d_{3}e^{cD_{e}}r_{\rm g}+d_{1}\bar{l}r_{\rm
g}\frac{L_{4}^{2}}{2}$ $\displaystyle+d_{1}r_{\rm g}\frac{L_{2}}{2}+r_{\rm
g}\frac{L_{1}}{2}+2d_{4}|P|+d_{2}r_{\rm g}L_{7}+d_{3}r_{\rm g}L_{8},$ (106)
$\displaystyle\beta_{2}=$ $\displaystyle d_{1}\frac{r_{\rm
g}^{2}}{2}+\frac{1024D_{e}^{2}}{D^{2}}r_{\rm
g}^{4}+8D_{e}^{2}\left(1+e^{-cD_{e}}\right)^{2}+\frac{1}{2\varepsilon}L_{5}\kappa$
$\displaystyle+d_{1}\frac{1}{\varepsilon}k_{n}L_{6}+\frac{64}{D}r_{\rm
g}^{2}\left(\bar{l}^{2}+\frac{\bar{l}^{3}}{3}\right)L_{4}^{2}+\frac{2}{d_{1}D}r_{\rm
g}^{2}L_{3}$
$\displaystyle+\frac{1}{2\epsilon}L_{5}\kappa+\frac{1}{2\epsilon}k_{n}L_{6}+2h_{n}+8k_{n}^{2}g^{2}.$
(107)

Taking the time derivative of Lyapunov functions in (97)-(100), using Agmon’s
inequality, Poincare’s inequality, Young’s inequality, and Assumptions 1, 2,
3, we obtain

$\displaystyle\dot{V}\leq$
$\displaystyle-\frac{D}{64}||\varpi_{xx}||^{2}-d_{1}\frac{D}{4}||\varpi_{x}||^{2}-d_{1}\frac{g}{16}||\varpi||^{2}-\frac{d_{2}}{2}|z(0,t)|^{2}$
$\displaystyle-\frac{d_{3}}{2}|z_{x}(0,t)|^{2}+d_{1}\frac{|\dot{l}(t)|^{2}}{2}||\varpi(t)||^{2}+\frac{8}{D}|\dot{l}(t)|^{2}|z(0,t)|^{2}$
$\displaystyle+d_{1}l(t)\left(\dot{l}(t)G_{1}(l(t))X(t)\int_{0}^{l(t)}\varpi(x,t)dx\right)$
$\displaystyle+d_{1}\dot{l}(t)\int_{0}^{l(t)}F(x,X(t))\varpi(x,t)dx$
$\displaystyle+\frac{1}{2\varepsilon}\int_{0}^{l(t)}\left(\phi(x-l(t))^{\top}f(X(t))\right)^{2}dx$
$\displaystyle+d_{1}\frac{1}{2\varepsilon}\int_{0}^{l(t)}\left(G_{2}(x-l(t))B\bar{h}(X)\right)^{2}dx$
$\displaystyle+\dot{l}(t)G_{1}(l(t)))X(t)\int_{0}^{l(t)}\left(l(t)+x\right)\varpi_{xx}(x,t)dx$
$\displaystyle+\frac{\big{|}\dot{l}(t)\big{|}}{2}F(l(t),X(t))^{2}$
$\displaystyle+\dot{l}(t)\int_{0}^{l(t)}F_{x}(x,X(t))\varpi_{x}(x,t)dx$
$\displaystyle+\frac{1}{2\epsilon}\int_{0}^{l(t)}\left(\phi(x-l(t))^{\top}f(X(t))\right)^{2}dx$
$\displaystyle+\frac{1}{2\epsilon}\int_{0}^{l(t)}\left(G_{3}(x-l(t))B\bar{h}(X)\right)^{2}dx$
$\displaystyle-d_{2}c\int_{0}^{D_{e}}e^{cx}z^{2}(x,t)dx-
d_{3}c\int_{0}^{D_{e}}e^{cx}z_{x}^{2}(x,t)dx$
$\displaystyle-2d_{2}\dot{l}(t)\int_{0}^{D_{e}}e^{cx}z(x,t)G_{4}(x-l(t))X(t)dx$
$\displaystyle-2d_{3}\dot{l}(t)\int_{0}^{D_{e}}e^{cx}z_{x}(x,t)G_{5}(x-l(t))X(t)dx$
$\displaystyle-d_{4}\frac{\lambda_{\rm
min}(Q)}{2}X^{\top}X+2d_{4}\big{|}P\big{|}\left(X^{\top}X\right)^{3/2}$
$\displaystyle+2\bar{h}^{\prime}(X)^{2}+2g^{2}\bar{h}(X)^{2}$ (108)

where

$\displaystyle G_{1}(l(t)):=\psi^{\prime}(-l(t))-q(0,l(t))C^{\top},$ (109)
$\displaystyle
G_{2}(x-l(t)):=-\phi^{\prime}(x-l(t))^{\top}-\frac{1}{D}\phi(x-l(t))^{\top},$
(110) $\displaystyle
G_{3}(x-l(t)):=-\phi^{\prime}(x-l(t))^{\top}-\frac{a}{D}\phi(x-l(t))^{\top},$
(111) $\displaystyle G_{4}(x-l(t)):=\psi^{\prime}(x-l(t))-q(x,l(t))C^{\top},$
(112) $\displaystyle
G_{5}(x-l(t)):=\psi^{\prime\prime}(x-l(t))-q_{x}(x,l(t))C^{\top},$ (113)

when we pick

$\displaystyle d_{1}\geq\frac{4a^{2}}{D^{2}},$ (114) $\displaystyle d_{2}\geq
d_{1}\left(g^{2}\frac{\bar{l}^{5}}{6D}+\frac{4(a^{2}+g^{2}\bar{l}^{2})+g}{g}\right)$
$\displaystyle\quad\quad+\frac{32(a^{2}+g^{2}\bar{l}^{2})}{D}+\frac{32\bar{l}^{3}}{3D}+d_{4}\frac{4\big{|}B^{T}P\big{|}^{2}}{\lambda_{\rm
min}(Q)},$ (115) $\displaystyle d_{3}\geq
d_{1}\frac{4\bar{l}^{2}}{g}+d_{1}\frac{4\bar{l}^{3}}{3g}+\frac{32\bar{l}^{5}}{3D}+\frac{32\bar{l}^{3}}{3D},$
(116) $\displaystyle d_{4}\leq\frac{D\lambda_{\rm
min}(Q)}{512\bar{l}\big{|}B^{T}P\big{|}^{2}}.$ (117)

With Assumption 1, there exist positive constants $L_{i}>0$ for $i=1,2,...,5$
such that the following inequalities hold:

$\displaystyle F(l(t),X(t))^{2}\leq L_{1}X^{\top}X,$ (118)
$\displaystyle\int_{0}^{l(t)}F(x,X(t))^{2}dx\leq L_{2}X^{\top}X,$ (119)
$\displaystyle\int_{0}^{l(t)}F_{x}(x,X(t))^{2}dx\leq L_{3}X^{\top}X,$ (120)
$\displaystyle\left|\psi^{\prime}(-l(t))-q(0,l(t))C^{\top}\right|\leq L_{4},$
(121) $\displaystyle\int_{0}^{l(t)}\phi(x-l(t))^{2}dx\leq L_{5},$ (122)
$\displaystyle\int_{0}^{l(t)}\left(G_{2}(x-l(t))B\right)^{2}dx\leq L_{6},$
(123) $\displaystyle\int_{0}^{l(t)}\left(G_{3}(x-l(t))B\right)^{2}dx\leq
L_{7}.$ (124)

In addition, $\dot{l}(t)$ can be rewritten as $\dot{l}(t)=r_{\rm g}e_{1}X$ by
(20). For the nonlinear terms, the upper bound for $f(X)$ is derived as

$\displaystyle f(X(t))\leq\kappa X^{\top}X.$ (125)

Next, we find the bound for (54). Since the constants defined in (10)-(11)
satisfy $\frac{a-gl_{\rm c}}{D}=\lambda_{+}K_{+}+\lambda_{-}K_{-}$, we can
show that (54) is equivalent to

$\displaystyle\bar{h}(X)=$ $\displaystyle
c_{\infty}K_{+}(1+\lambda_{+}z_{2}-e^{\lambda_{+}z_{2}})$
$\displaystyle+c_{\infty}K_{-}(1+\lambda_{-}z_{2}-e^{\lambda_{-}z_{2}}).$
(126)

Since we know that $e^{x}\leq 1+x+x^{2}$ for $x<1.79$, we can show the bound
of (126) for local region as

$\displaystyle\left|\bar{h}(X)\right|$ $\displaystyle\leq
c_{\infty}(\left|K_{+}\right|(\lambda_{+}z_{2})^{2}+\left|K_{-}\right|(\lambda_{-}z_{2})^{2}).$
(127)

Now, we can pick
$k_{n}=\max\\{c_{\infty}K_{+}\lambda_{+}^{2},c_{\infty}K_{-}\lambda_{-}^{2}\\}$.
Thus,

$\displaystyle\left|\bar{h}(X)\right|$ $\displaystyle\leq
2\left|k_{n}\right|z_{2}^{2}\leq 2k_{n}X^{\top}X.$ (128)

Similarly, the upper bound for $\bar{h}^{\prime}(X)$ term is obtained for
$x<1.79$ as

$\displaystyle\bar{h}^{\prime}(X)^{2}$ $\displaystyle\leq
h_{n}(X^{\top}X)^{2}$ (129)

where $h_{n}=\max\left\\{2c_{\infty}r_{\rm
g}^{2}K_{+}^{2}\lambda_{+}^{4},2c_{\infty}r_{\rm
g}^{2}K_{-}^{2}\lambda_{-}^{4}\right\\}$.

By using Cauchy-Schwarz, Agmon’s, Young’s inequalities, and the following
inequality

$\displaystyle||z_{x}(x,t)||^{2}$
$\displaystyle\leq\int_{0}^{D_{e}}e^{cx}z_{x}(x,t)^{2}dx\leq
e^{cD_{e}}||z_{x}(x,t)^{2}||,$ (130)

we obtain

$\displaystyle\dot{V}\leq$ $\displaystyle-
d_{1}\frac{D}{8}V_{2}-d_{1}\frac{g}{16}V_{1}-d_{4}\frac{\lambda_{\rm
min}(Q)}{2\lambda_{\min}(P)}V_{5}-d_{2}cV_{3}$ $\displaystyle-
d_{3}cV_{4}+\left(d_{1}\frac{r_{\rm
g}\left(\bar{l}+1\right)}{2\lambda_{\min}(P)}+\frac{e^{cD_{e}}r_{\rm
g}\left(d_{2}+d_{3}\right)}{\lambda_{\min}(P)}\right)V^{3/2}$
$\displaystyle+\left(\frac{d_{1}r_{\rm
g}L_{2}}{2\lambda_{\min}(P)}+\frac{r_{\rm
g}L_{1}}{2\lambda_{\min}(P)}++\frac{d_{1}\bar{l}r_{\rm
g}L_{4}^{2}}{2\lambda_{\min}(P)}\right)V^{3/2}$
$\displaystyle+\frac{1}{\lambda_{\min}(P)}\left(2d_{4}|P|+d_{2}r_{\rm
g}L_{7}+d_{3}r_{\rm g}L_{8}\right)V^{3/2}$
$\displaystyle+\frac{1}{\lambda_{\min}(P)}\left(d_{1}\frac{r_{\rm
g}^{2}}{2}+\frac{1024D_{e}^{2}}{D^{2}}r_{\rm
g}^{4}+8cD_{e}^{4}+\frac{L_{5}}{2\varepsilon}\kappa\right)V^{2}$
$\displaystyle+\frac{1}{\lambda_{\min}(P)}\left(d_{1}\frac{L_{6}+L_{7}}{\varepsilon}k_{n}+\frac{64}{D}r_{\rm
g}^{2}\left(\bar{l}^{2}+\frac{\bar{l}^{3}}{3}\right)L_{4}^{2}\right)V^{2}$
$\displaystyle+\frac{1}{\lambda_{min}(P)}\left(\frac{2r_{\rm
g}^{2}L_{3}}{d_{1}D}+\frac{L_{5}}{2\epsilon}\kappa+\frac{L_{6}}{2\epsilon}k_{n}+2h_{n}\right)V^{2}$
$\displaystyle+\frac{1}{\lambda_{\min}(P)}8k_{n}^{2}g^{2}V^{2}.$ (131)

The inequality above can be written as

$\displaystyle\dot{V}$ $\displaystyle\leq-\alpha
V+\beta_{1}V^{3/2}+\beta_{2}V^{2},$ (132)

where $\alpha$, $\beta_{1}$ and $\beta_{2}$ are defined in (105)-(107), which
completes the proof of Lemma 1.

###  4.2 Guaranteeing the conditions for all time

Next, we prove the following lemmas to conclude Theorem 1 by ensuring the
local stability of the closed-loop system.

###### Lemma 2

Under Assumption 3, there exists a positive constant $K_{1}>0$ such that in
$\Omega_{1}:=\\{(\varpi,z,X)\in H_{1}\times
H_{1}\times\mathbb{R}^{2}|V(t)<K_{1}\\}$, Assumptions 1 and 2 are satisfied.

The definition of $X(t)$ in (25) is equivalent to

$\displaystyle X(t)=\left[\begin{array}[]{cc}\frac{\dot{l}(t)}{r_{\rm
g}}&l(t)-l_{\rm s}\end{array}\right]^{\top}.$ (134)

Thus, for some $r>0$, if $|X|<r$ then both the following two inequalities
hold:

$\displaystyle\left|\frac{\dot{l}(t)}{r_{\rm g}}\right|<r,\ \
\big{|}l(t)-l_{\rm s}\big{|}<r.$ (135)

The first inequality ensures that if $r<\frac{\bar{v}}{r_{\rm g}}$ then
Assumption 1 holds. The second inequality yields $-r+l_{\rm s}<l(t)<r+l_{\rm
s}$, and thus if both $r<l_{\rm s}$ and $r<\bar{l}-l_{\rm s}$ hold, then
Assumption 2 holds. Hence, $r>0$ is chosen as

$\displaystyle r=\min\left\\{\frac{\bar{v}}{r_{\rm g}},l_{\rm
s},\bar{l}-l_{\rm s}\right\\}.$ (136)

Since the norm of $X$ is bounded as

$\displaystyle|X|^{2}\leq\frac{1}{\lambda_{\rm
min}(P)}X^{\top}PX\leq\frac{d_{4}}{\lambda_{\rm min}(P)}V.$ (137)

Then, by setting $K_{1}=\frac{\lambda_{\rm min}(P)}{d_{4}}r^{2}$, if
$V(t)<K_{1}$ holds then $|X|<r$ and thus Assumptions 1 and 2 are satisfied, by
which we can conclude Lemma 2.

###### Lemma 3

Under Assumption 3, there exists a positive constant $K_{2}>0$ such that if
$V(0)<K_{2}$ then the Assumption 1 and 2 are satisfied and the following norm
estimate holds:

$\displaystyle V(t)\leq V(0)e^{-\frac{\alpha}{2}t}.$ (138)

For a positive constant $K_{2}>0$, let $\Omega:=\\{(\varpi,z,X)\in H_{1}\times
H_{1}\times\mathbb{R}^{2}|V(t)<K_{2}\\}$. If $K_{2}\leq K_{1}$ then
$\Omega\subset\Omega_{1}$, and thus the assumptions are satisfied in $\Omega$.
Moreover, due to Lemma 1, the norm estimate (104) holds. Hence, by setting

$\displaystyle
K_{2}\leq\frac{-\beta_{2}^{2}\sqrt{\frac{\beta_{1}^{2}(2\alpha\beta_{2}+\beta_{1}^{2})}{\beta_{2}^{4}}}+\alpha\beta_{2}+\beta_{1}^{2}}{2\beta_{2}^{2}},$
(139)

we can see that applying $V(t)<K_{2}$ to (104) leads to

$\displaystyle\dot{V}\leq-\frac{\alpha}{2}V,$ (140)

by which the norm estimate (138) is deduced. Since (138) is a monotonically
decreasing function in time, by setting $K_{3}=\min\\{K_{1},K_{2}\\}$, the
region $\Omega$ is shown to be an invariant set. Thus, if $V(0)<K_{3}$, then
$V(t)<K_{3}$ for all $t\geq 0$, and one can conclude with Lemma 3. Finally, we
show that $(w,z,X)$ system is locally exponentially stable. Taking the square
of (89) and applying Young’s and Cauchy-Schwarz inequalities, we obtain

$\displaystyle||\varpi(\cdot,t)||_{H_{1}}^{2}\leq
3||w(\cdot,t)||_{H_{1}}^{2}+M_{1}||z(\cdot,t)||_{H_{1}}^{2}+M_{2}|X|^{2},$
(141) $\displaystyle||w(\cdot,t)||_{H_{1}}^{2}\leq
3||\varpi(\cdot,t)||_{H_{1}}^{2}+M_{1}||z(\cdot,t)||_{H_{1}}^{2}+M_{2}|X|^{2},$
(142)

where $M_{1}=4D_{e}^{3}\bar{l}+4D_{e}\bar{l}^{3}$ and $M_{2}=12k_{n}^{2}$.
Consider the norm
$\Pi(x,t)=||w(x,t)||_{H_{1}}^{2}+||z(x,t)||_{H_{1}}^{2}+|X(t)|^{2}$. By using
(130), we can show that

$\displaystyle||z(.,t)||_{H_{1}}^{2}\leq V_{3}+V_{4}\leq
e^{cD_{e}}||z(.,t)||_{H_{1}}^{2},$ (143)

and

$\displaystyle M_{3}||w(.,t)||_{H_{1}}^{2}\leq 2V_{1}+2V_{2}\leq
M_{4}||w(.,t)||_{H_{1}}^{2},$ (144)

where $M_{3}=\max\\{\bar{l}^{2},1\\}$ and
$M_{4}=\max\\{\frac{1}{\bar{l}^{2}},1\\}$. Then, we have

$\displaystyle\Pi(t)\leq$
$\displaystyle\left(1+M_{1}\right)||z_{x}(x,t)|||_{H_{1}}^{2}+3||\varpi(x,t)||_{H_{1}}^{2}$
$\displaystyle+(1+M_{2})|X(t)|^{2}$ $\displaystyle\leq$
$\displaystyle(2+2M_{1})(V_{3}+V_{4})+12(V_{1}+V_{2})+\frac{1+M_{2}}{\lambda_{\min}(P)}V_{5}$
(145)

By using (96), we have

$\displaystyle V\leq$
$\displaystyle\left(\max\\{d_{2},d_{3}\\}e^{cD_{e}}+\frac{M_{1}M_{3}}{2}\max\\{d_{1},1\\}\right)||z(.,t)||_{H_{1}}^{2}$
$\displaystyle+\frac{3M_{3}}{2}\max\\{d_{1},1\\}||w(.,t)||_{H_{1}}^{2}$
$\displaystyle+\left(\frac{d_{4}}{\lambda_{\min}(P)}+\frac{M_{2}M_{3}}{2}\max\\{d_{1},1\\}\right)|X(t)|^{2}$
$\displaystyle\leq$
$\displaystyle\Sigma_{1}||z(.,t)||_{H_{1}}^{2}+\Sigma_{2}||w(.,t)||_{H_{1}}^{2}+\Sigma_{3}|X(t)|^{2}$
(146)

where
$\Sigma_{1}=\max\\{d_{2},d_{3}\\}e^{cD_{e}}+\frac{M_{1}M_{3}}{2}\max\\{d_{1},1\\}$,
$\Sigma_{2}=\frac{3M_{3}}{2}\max\\{d_{1},1\\}$ and
$\Sigma_{3}=\frac{d_{4}}{\lambda_{\min}(P)}+\frac{M_{2}M_{3}}{2}\max\\{d_{1},1\\}$.
Therefore, (145) and (146) leads us the following norm equivalence

$\displaystyle\underline{$\delta$}V(t)\leq\Pi(t)\leq\bar{\delta}V(t)$ (147)

where

$\displaystyle\underline{$\delta$}$
$\displaystyle=\frac{1}{\max\left\\{\Sigma_{1},\Sigma_{2},\Sigma_{3}\right\\}},$
(148) $\displaystyle\bar{\delta}$
$\displaystyle=\max\left\\{(2+2M_{1}),12,\frac{1+M_{2}}{\lambda_{\min}(P)}\right\\}.$
(149)

With the inequalities (138) and (147), we show that

$\displaystyle\Pi(t)\leq\frac{\bar{\delta}}{\underline{$\delta$}}\Pi(0)e^{-\frac{\alpha}{2}},$
(150)

which ensures the exponential stability of $(w,z,X)$-system. Finally, owing to
Lemma 3 and (150), and using the similar norm equivalent estimate in the
$H_{1}$-norm between the target system $(w,z,X)$ and the closed-loop system
$(u,v,X)$, the local stability of the closed-loop system is proved, which
completes the proof of Theorem 1.

##  5 Conclusion

This paper proposes a novel feedback control law with input delay compensation
for a PDE with a moving boundary governed by an ODE. The input delay is
represented by a transport PDE dynamics, and the controller is designed by PDE
backstepping for stabilizing the PDE-ODE system with compensating the input
delay. For further research, the observer design can be studied under delayed
outputs so that the unknown state variables can be estimated through
compensating the delayed measured output.

## References

  * Bradbury and Carter (2011) Bradbury, E.J. and Carter, L.M. (2011).  Manipulating the glial scar: chondroitinase abc as a therapy for spinal cord injury.  _Brain research bulletin_ , 84(4-5), 306–316.
  * Buisson-Fenet et al. (2018) Buisson-Fenet, M., Koga, S., and Krstic, M. (2018).  Control of piston position in inviscid gas by bilateral boundary actuation.  In _2018 IEEE Conference on Decision and Control (CDC)_ , 5622–5627. IEEE.
  * Clark et al. (2009) Clark, B.D., Goldberg, E.M., and Rudy, B. (2009).  Electrogenic tuning of the axon initial segment.  _The Neuroscientist_ , 15(6), 651–668.
  * Demir et al. (2021) Demir, C., Koga, S., and Krstic, M. (2021).  Neuron growth control by pde backstepping: Axon length regulation by tubulin flux actuation in soma.  In _2021 60th IEEE Conference on Decision and Control (CDC)_ , 649–654.  10.1109/CDC45484.2021.9683188.
  * Demir et al. (2022) Demir, C., Koga, S., and Krstic, M. (2022).  Neuron growth output-feedback control by pde backstepping.  _arXiv preprint arXiv:2203.16733_.
  * Desai and Mitchison (1997) Desai, A. and Mitchison, T.J. (1997).  Microtubule polymerization dynamics.  _Annual review of cell and developmental biology_ , 13(1), 83–117.
  * Diehl et al. (2014) Diehl, S., Henningsson, E., Heyden, A., and Perna, S. (2014).  A one-dimensional moving-boundary model for tubulin-driven axonal growth.  _Journal of theoretical biology_ , 358, 194–207.
  * Frantz et al. (2010) Frantz, C., Stewart, K.M., and Weaver, V.M. (2010).  The extracellular matrix at a glance.  _Journal of cell science_ , 123(24), 4195–4200.
  * García-Grajales et al. (2017) García-Grajales, J.A., Jérusalem, A., and Goriely, A. (2017).  Continuum mechanical modeling of axonal growth.  _Computer Methods in Applied Mechanics and Engineering_ , 314, 147–163.
  * Julien (1999) Julien, J.P. (1999).  Neurofilament functions in health and disease.  _Current opinion in neurobiology_ , 9(5), 554–560.
  * Kandel et al. (2000) Kandel, E.R., Schwartz, J.H., Jessell, T.M., Siegelbaum, S., Hudspeth, A.J., and Mack, S. (2000).  _Principles of neural science_ , volume 4.  McGraw-hill New York.
  * Karimi-Abdolrezaee et al. (2010) Karimi-Abdolrezaee, S., Eftekharpour, E., Wang, J., Schut, D., and Fehlings, M.G. (2010).  Synergistic effects of transplanted adult neural stem/progenitor cells, chondroitinase, and growth factors promote functional repair and plasticity of the chronically injured spinal cord.  _Journal of Neuroscience_ , 30(5), 1657–1676.
  * Koga et al. (2020) Koga, S., Bresch-Pietri, D., and Krstic, M. (2020).  Delay compensated control of the stefan problem and robustness to delay mismatch.  _International Journal of Robust and Nonlinear Control_ , 30(6), 2304–2334.
  * Koga et al. (2018) Koga, S., Diagne, M., and Krstic, M. (2018).  Control and state estimation of the one-phase stefan problem via backstepping design.  _IEEE Transactions on Automatic Control_ , 64(2), 510–525.
  * Koga and Krstic (2020) Koga, S. and Krstic, M. (2020).  _Materials Phase Change PDE Control and Estimation: From Additive Manufacturing to Polar Ice_.  Springer Nature.
  * Krstic (2009a) Krstic, M. (2009a).  Control of an unstable reaction–diffusion pde with long input delay.  _Systems & Control Letters_, 58(10-11), 773–782.
  * Krstic (2009b) Krstic, M. (2009b).  Delay compensation for nonlinear, adaptive, and pde systems.
  * Krstic and Smyshlyaev (2008) Krstic, M. and Smyshlyaev, A. (2008).  _Boundary control of PDEs: A course on backstepping designs_.  SIAM.
  * Lee et al. (2010) Lee, H., McKeon, R.J., and Bellamkonda, R.V. (2010).  Sustained delivery of thermostabilized chabc enhances axonal sprouting and functional recovery after spinal cord injury.  _Proceedings of the National Academy of Sciences_ , 107(8), 3340–3345.
  * Lemons et al. (1999) Lemons, M.L., Howland, D.R., and Anderson, D.K. (1999).  Chondroitin sulfate proteoglycan immunoreactivity increases following spinal cord injury and transplantation.  _Experimental neurology_ , 160(1), 51–65.
  * Liu et al. (1997) Liu, X.Z., Xu, X.M., Hu, R., Du, C., Zhang, S.X., McDonald, J.W., Dong, H.X., Wu, Y.J., Fan, G.S., Jacquin, M.F., et al. (1997).  Neuronal and glial apoptosis after traumatic spinal cord injury.  _Journal of Neuroscience_ , 17(14), 5395–5406.
  * Maccioni et al. (2001) Maccioni, R.B., Muñoz, J.P., and Barbeito, L. (2001).  The molecular bases of alzheimer’s disease and other neurodegenerative disorders.  _Archives of medical research_ , 32(5), 367–381.
  * McLean and Graham (2006) McLean, D.R. and Graham, B.P. (2006).  Stability in a mathematical model of neurite elongation.  _Mathematical medicine and biology: a journal of the IMA_ , 23(2), 101–117.
  * McLean et al. (2004) McLean, D.R., van Ooyen, A., and Graham, B.P. (2004).  Continuum model for tubulin-driven neurite elongation.  _Neurocomputing_ , 58, 511–516.
  * Smyshlyaev and Krstic (2004) Smyshlyaev, A. and Krstic, M. (2004).  Closed-form boundary state feedbacks for a class of 1-d partial integro-differential equations.  _IEEE Transactions on Automatic Control_ , 49(12), 2185–2202.
  * Susto and Krstic (2010) Susto, G.A. and Krstic, M. (2010).  Control of pde–ode cascades with neumann interconnections.  _Journal of the Franklin Institute_ , 347(1), 284–314.
  * Tang and Xie (2011) Tang, S. and Xie, C. (2011).  State and output feedback boundary control for a coupled pde–ode system.  _Systems & Control Letters_, 60(8), 540–545.
  * Tessier-Lavigne and Goodman (1996) Tessier-Lavigne, M. and Goodman, C.S. (1996).  The molecular biology of axon guidance.  _Science_ , 274(5290), 1123–1133.
  * Van Veen and Van Pelt (1994) Van Veen, M.P. and Van Pelt, J. (1994).  Neuritic growth rate described by modeling microtubule dynamics.  _Bulletin of mathematical biology_ , 56(2), 249–273.
  * Yamada et al. (1970) Yamada, K.M., Spooner, B.S., and Wessells, N.K. (1970).  Axon growth: roles of microfilaments and microtubules.  _Proceedings of the National Academy of Sciences_ , 66(4), 1206–1212.
  * Yu et al. (2020) Yu, H., Diagne, M., Zhang, L., and Krstic, M. (2020).  Bilateral boundary control of moving shockwave in lwr model of congested traffic.  _IEEE Transactions on Automatic Control_.


#  Spatial organization of slit-confined melts of ring polymers with non-
conserved topology: A lattice Monte Carlo study

Mattia Alberto Ubertini<EMAIL_ADDRESS>Scuola Internazionale Superiore di
Studi Avanzati (SISSA), Via Bonomea 265, 34136 Trieste, Italy  Angelo Rosa
<EMAIL_ADDRESS>Scuola Internazionale Superiore di Studi Avanzati (SISSA), Via
Bonomea 265, 34136 Trieste, Italy

###### Abstract

We present Monte Carlo computer simulations for melts of semiflexible randomly
knotted and randomly concatenated ring polymers on the fcc lattice and in slit
confinement. Through systematic variation of the slit width at fixed melt
density, we first explore the influence of confinement on single-chain
conformations and inter-chain interactions. We demonstrate that confinement
makes chains globally larger and more elongated, while enhancing both contacts
and knottedness propensities. As for multi-chain properties, we show that
ring-ring contacts decrease with the confinement, yet neighbouring rings are
more overlapped as confinement grows. These aspects are reflected on the
decrease of the links formation between pairs of rings. The results suggest
that confinement can be used to fine-tune the mechanical properties of the
polymer network. In particular, confinement biases the synthesis of networks
that are softer to mechanical stress. Finally, in connection with a previous
study of us and recent simulations on two-dimensional polymer melts, our
findings suggest that entanglements in polymer melts arise from pairwise ring-
ring links alone.

##  I 1\. Introduction

Recent years have witnessed a growing interest in the design of so called
smart materials, such polycatenanes and polyrotaxanes Wu et al. (2017); Hart
et al. (2021), whose microscopic components are constituted by ring polymers
interlocked to each other by topological links that can be artificially
synthesised following precise chemical routes. Interestingly similar devices
can be also prepared by employing biological components, mainly DNA plasmid
rings Krajina et al. (2018) which interlock to each other through the action
of the enzyme topoisomerase-II (TopoII) and form a molecular state termed
Olympic hydrogel which has been first theorized by de Gennes in 1997 Raphaël
et al. (1997). Remarkably, similar molecules can be also found in Nature: a
classical example is the kinetoplast DNA Chen et al. (1995) present in the
mitochondria of certain Trypanosoma parasites.

Similarly to covalent bonds stabilizing the shape of a molecule, topological
links remain stable at room temperature which guarantees the corresponding
molecule to maintain a relatively well characterized spatial conformation. On
the other hand, since the single ring constituents are not rigid objects but
they fluctuate Chiarantoni and Micheletti (2022) as ordinary polymers
typically do Doi and Edwards (1986); Rubinstein and Colby (2003), these
molecules display unusual mechanical properties under stress and tunable
viscoelasticity that can be exploited in a wide number of practical
applications (molecular machines and drug delivery Rauscher et al. (2020a, b),
to name a few), thus justifying the adjective “smart” employed for these
materials.

The preparation of topological materials with well designed properties is a
delicate balance between many parameters: indeed, several numerical studies
Lang et al. (2012, 2014); Ubertini and Rosa (2021, 2023) have characterized
the topological state of systems made up of randomly concatenated and knotted
polymer rings, and have shown that the resulting networks can be controlled
using experimentally tunable parameters such as the length of the polymer
chain, the density of the polymer solution, and the bending stiffness of the
polymer fiber. So far though, geometric confinement as a way to drive the
synthesis of concatenated ring networks has received considerable less
attention. Yet, recent experiments Soh and Doyle (2021) performed on
kinetoplast DNA Chen et al. (1995) at varying degree of slit confinement have
foreseen the possibility of exploiting geometric constraints to bias the
synthesis of a DNA-based network, similarly to the one discussed in Ref.
Krajina et al. (2018).

In this work, we explore how geometric constraints, under the form of slit
confinement, can affect the structural properties of systems of strand-
crossing rings. To this purpose, we perform extensive dynamical simulations of
highly entangled systems of randomly concatenated and knotted rings employing
the kinetic Monte Carlo algorithm introduced by us Ubertini and Rosa (2021)
for studying these systems at bulk conditions. Varying the degree of
confinement, we quantify its influence on the metric properties of the rings,
which present interesting non-monotonous behaviour, as well as topological
ones, in particular knotting probability is highly enhanced by reducing the
height of the slit, while the linking between the rings is diminished. These
findings suggest that geometric confinement can be used as a powerful tool to
control the topology of the resulting networks and their elastic properties.

The paper is structured as the following. In Section 1 we present and discuss
the Monte Carlo lattice polymer model, we introduce the notation and we
explain how to detect and compute topological invariants for the
characterization of knots and links in the system. In Sec. 2 we present the
main results of our work, while in Sec. 3 we provide some discussion and
conclusions regarding the role of slit confinement in shaping both single-
chain and inter-chain properties of the resulting polymer networks. Additional
figures have been included in the Supporting Information (SI) file.

##  II 2\. Model and methods

###  II.1 2.1. Polymer model

We consider polymer melts made of $M$ randomly concatenated and randomly
knotted ring polymers of $N=320$ monomers each on the fcc lattice; the fcc
unit step $a$ is taken as our unit length. The simulations are based on the
kinetic Monte Carlo (kMC) algorithm introduced by us in Ubertini and Rosa
(2021). Since then, the algorithm has been variously applied to study melts of
non-concatenated and unknotted rings Ubertini et al. (2022) and the connection
between entanglements and physical links in semiflexible chain melts Ubertini
and Rosa (2023). In this article we limit ourselves to summarizing the
essential details of the numerical protocol, while referring the reader to our
past works for more details.

${\rm H}/a$ |  $\hat{{\rm H}}$ |  $M$ |  $\langle b\rangle/a$
---|---|---|---
$2.12$ | $0.30$ | $420$ | $0.656$
$3.53$ | $0.50$ | $422$ | $0.658$
$4.95$ | $0.70$ | $420$ | $0.659$
$6.36$ | $0.90$ | $427$ | $0.659$
$7.78$ | $1.10$ | $420$ | $0.660$
$10.61$ | $1.51$ | $420$ | $0.660$
$13.43$ | $1.91$ | $422$ | $0.660$
$17.68$ | $2.51$ | $430$ | $0.660$
$20.51$ | $2.91$ | $433$ | $0.660$
${\rm bulk}$ | – | $420$ | $0.663$
Table 1:  Values of physical parameters for the ring polymer melts
investigated in this paper. $a$ is the unit distance of the fcc lattice and
the monomer number per fcc lattice site is equal to $\frac{5}{4}=1.25$, see
text and Refs. Ubertini and Rosa (2021); Ubertini et al. (2022); Ubertini and
Rosa (2023) for details. (i) ${\rm H}$, height of the slit. (ii) $\hat{\rm
H}=\frac{{\rm H}}{\sqrt{\langle R^{2}_{g}\rangle_{\rm bulk}}}$, ratio between
the height of the slit and the root-mean-square gyration radius of rings in
bulk (i.e., no confinement) conditions. (iii) $M$, total number of simulated
chains in the melt. (iv) $\langle b\rangle$, mean bond length Mea .

Essentially, the polymer model takes into account: (i) chain connectivity,
(ii) bending stiffness, (iii) excluded volume, (iv) topological rearrangement
of polymer chains. Finally, and for the first time, in this work we consider
(v) slit confinement in the model. For the implementation of chain dynamics,
the following combination of MC moves – that automatically take into account
excluded volume interactions – are used:

  * (a)

Topology-preserving moves (termed Rouse-like and reptation-like, see Ubertini
and Rosa (2021)) that automatically enforce excluded volume interactions. By
construction, these moves enable two (and no more than two) consecutive bonded
monomers along each single chain to occupy the same lattice site: by allowing
to store contour length along the polymer filament, this numerical “trick”
makes the chains locally elastic and facilitates global chain equilibration.
Because of that, the bond length is a fluctuating quantity with mean value
$=\langle b\rangle$: in particular, the latter is insensitive to confinement
(the measured values for $\langle b\rangle$ are reported in Table 1). In this
way, the mean polymer contour length is $L=N\langle b\rangle$ and, similarly,
the mean contour length of a sub-chain of $n$ monomers is $\ell=n\langle
b\rangle$.

  * (b)

Topology-changing moves Ubertini and Rosa (2021) that induce random strand
crossings between nearby polymer filaments at a tunable rate: we set this rate
to $10^{4}$ kMC elementary steps, consistent with our previous works Ubertini
and Rosa (2021); Ubertini et al. (2022); Ubertini and Rosa (2023). Strand-
crossings between filaments of the same ring can result in the creation or
destruction of knots, while inter-ring crossings may cause either catenation
or decatenation. The model has been shown to exhibit dynamical behavior
consistent with the experiments Krajina et al. (2018), specifically dynamic
“fluidization” of the rings due to topological violations through strand-
crossings. Thus, by performing simulations of strand-crossing rings, we sample
the ensemble of the network structures formed by randomly concatenated and
knotted rings at the given density and in slit confinement (see below for
details).

Then, bending stiffness is modelled in terms of the Hamiltonian (in Boltzmann
units, $\kappa_{B}T$):

$\frac{\mathcal{H}}{\kappa_{B}T}=-\kappa_{\rm bend}\\!\sum_{i=1}^{N\\!\langle
b\rangle/a}\\!\cos\theta_{i}\,,$

where $\kappa_{\rm bend}=2$ is the bending stiffness and $\theta_{i}$ is the
angle between consecutive bonds along the chain, with periodic conditions –
due to ring geometry – assumed for all the chains. By fixing the monomer
number per fcc lattice site equal to $\frac{5}{4}=1.25$ Ubertini and Rosa
(2021); Ubertini et al. (2022); Ubertini and Rosa (2023), the chosen bending
stiffness corresponds to the chain Kuhn segment $\ell_{K}/a=3.4$ Ubertini et
al. (2022) which is high enough to guarantee that distinct polymers are in an
effective highly entangled state.

Finally, ring polymers are subject to slit confinement. This particular form
of constraint is imposed by forcing the chains to move on the fcc lattice,
with periodic boundary conditions on the $xy$-plane and hard boundaries in the
$z$-direction placed in $z=0$ and $z={\rm H}$. We vary the height of the box
${\rm H}$ to study different confinement regimes, while adjusting the lateral
box sides $L_{x}=L_{y}$ to keep density constant. The degree of confinement is
quantified by the ratio $\hat{\rm H}={\rm H}/\sqrt{\langle
R^{2}_{g}\rangle_{\rm bulk}}$, expressing the ratio between the height, or
width, of the slit $\rm H$ and the root-mean-square gyration radius (see
definition (3)), $\sqrt{\langle R^{2}_{g}\rangle_{\rm bulk}}/a=\sqrt{49.66}$
Ubertini et al. (2022), of rings in bulk conditions. We investigate system’s
behavior from highly confined ($\hat{\rm H}\simeq 0.30$) to mildly confined
($\hat{\rm H}\simeq 2.91$) regimes and systematically compare the results with
the corresponding values in bulk. Wherever appropriate, we have also compared
the systems here with melts of unknotted and non-concatenated rings in bulk
Ubertini et al. (2022). We simulate $M\simeq 420$ chains, comprising a total
of $N\times M\simeq 134400$ monomers, with $M$ slightly adjusted to maintain a
constant density (see Table 1 for specific numbers).

To assess meaningful chain statistics and as in our other works Ubertini and
Rosa (2021, 2023) on similar polymer systems, we run simulations long enough
in order to get properly equilibrated melts. This is visualized in Figure S1
in SI, that shows plots of the monomer time mean-square displacement in the
frame of the centre of mass of the corresponding chain (the so called $g_{2}$
Kremer and Grest (1990)) as a function of the MC simulation time $\tau_{\rm
MC}$. As known, provided long-enough simulations are available, $g_{2}$
displays a plateau that is indicative of the equilibration of the system. All
our systems display corresponding plateaus, that demonstrates that
equilibration has been reached for all the cases considered. Accordingly, the
time scale to reach the corresponding plateau corresponds to the portion of
the trajectory that has been discarded from the computation of the relative
observables.

###  II.2 2.2. Detection of knots and links

In order to characterize the topological states of the rings in the melt, we
follow closely the pipeline recently developed by us Ubertini and Rosa (2023).
Specifically, we employ a numerical algorithm which “shrinks” or simplifies
each ring to its “primitive” shape, i.e. without violating topological
constraints: in this way we detect knots and links at any order, i.e. pairwise
links as well as three-chain links like the Borromean ring configuration
$6^{3}_{2}$ (see Sec. 2.3). The algorithm is able to return the irreducible
knotted or linked structure which we further characterize by computing their
topological invariants. For knots, in particular, we compute the corresponding
Jones polynomial Jones (1985) using the Python package Topoly Dabrowski-
Tumanski et al. (2021). Instead, for two-body links we compute the Gauss
linking number (GLN):

${\rm
GLN}\equiv\frac{1}{4\pi}\oint_{{\mathcal{C}}_{1}}\oint_{{\mathcal{C}}_{2}}\frac{(\vec{r}_{2}-\vec{r}_{1})\cdot(d{\vec{r}}_{2}\wedge
d{\vec{r}}_{1})}{|\vec{r}_{2}-\vec{r}_{1}|^{3}}\,,$ (1)

which gives the number of times two closed loops $\mathcal{C}_{1}$ and
$\mathcal{C}_{2}$, parametrized respectively by coordinates $\vec{r}_{1}$ and
$\vec{r}_{2}$, wind around each other. While unconcatenated rings have
$\rm{GLN}=0$, it is known that there exists concatenated pairs with
$\rm{GLN}=0$ (for instance, the so called Whitehead link configuration
$5^{2}_{1}$). In these “pathological” cases, the ones detected via our
shrinking algorithm have been successively identified by computing the Jones
polynomial using Topoly again. We compute the Jones polynomials also for
three-body irreducible links (for instance, Borromean rings) where a pairwise
topological invariant such as the ${\rm GLN}$ fails (see Sec. 3.2.2).

###  II.3 2.3. Notation

As for rings’ metric properties, for some observables $\mathcal{O}$ which can
be expressed as a function of monomers’ coordinates we study separately the
contributions $\mathcal{O}_{\perp}$ and $\mathcal{O}_{\parallel}$,
respectively perpendicular (or, transverse) and parallel to the plane of the
slit (which, by construction (see Sec. 2.1), coincides with the $xy$-plane).

As for rings’ topological properties, in referring to a given knot or link we
employ the conventional notation illustrated in the book by Rolfsen Rolfsen
(2003). Namely, a knot or a link is defined by the symbol $K_{i}^{p}$ where:
$K$ represents the number of irreducible crossings of the knot (or the link),
$p$ is the number of rings which takes part in the topological structure
(e.g., $p=2$ for two-chain links) and $i$ is an enumerative index assigned to
distinguish topologically non-equivalent structures having the same $K$ and
$p$.

##  III 3\. Results

###  III.1 3.1. Single-chain properties

####  III.1.1 3.1.1. Rings’ size and shape

First, we characterize the impact of slit confinement on the size and shape of
the rings. To this purpose, for each ring of the system we compute the
$3\times 3$ symmetric gyration tensor $Q_{\alpha\beta}=Q_{\beta\alpha}$
($\alpha,\beta=x,y,z$) defined as:

$Q_{\alpha\beta}=\frac{1}{N}\sum_{m=1}^{N}\left(r_{m,\alpha}-r_{{\rm
CM},\alpha}\right)\left(r_{m,\beta}-r_{{\rm CM},\beta}\right)\,,$ (2)

where $r_{m,\alpha}$ is the $\alpha$-th Cartesian component of the spatial
position $\vec{r}_{m}$ of monomer $m$ and $\vec{r}_{\rm
CM}\equiv\frac{1}{N}\sum_{m=1}^{N}\vec{r}_{m}$ is the center of mass of the
chain. The mean values of the eigenvectors of $Q$ ordered in descending order,
$\langle\lambda^{2}_{1}\rangle\geq\langle\lambda^{2}_{2}\rangle\geq\langle\lambda^{2}_{3}\rangle$,
quantify the mean spatial elongations of the polymers on the corresponding
principal axes, while the mean value of the trace of $Q$, $\langle{\rm
tr}Q\rangle=\sum_{\alpha=1}^{3}\langle\lambda_{\alpha}^{2}\rangle$, is equal
to the mean-square gyration radius or size,

$\langle
R_{g}^{2}\rangle\equiv\frac{1}{N}\sum_{m=1}^{N}\langle\left({\vec{r}}_{m}-{\vec{r}}_{\rm
CM}\right)^{2}\rangle=\langle{\rm
tr}Q\rangle=\sum_{\alpha=1}^{3}\langle\lambda_{\alpha}^{2}\rangle\,,$ (3)

of the chain.

Figure 1:  (a) Ring mean-square gyration radius ($\langle R_{g}^{2}\rangle$)
with its parallel ($\langle R_{g,\parallel}^{2}\rangle$) and transverse
($R_{g,\perp}^{2}$) components as a function of the degree of confinement
$\hat{\rm H}$ (see Sec. 2.1 for definition). The dashed lines are for the
values of the bulk system (i.e., no confinement). Error bars are smaller than
the symbols size. (b, c, d) Scaling plots for, respectively, distribution
functions of the ring gyration radius ($P(R_{g}/\sqrt{\langle
R_{g}^{2}\rangle})$) and of its parallel ($P(R_{g,\parallel}/\sqrt{\langle
R_{g,\parallel}^{2}\rangle})$) and transverse ($P(R_{g,\perp}/\sqrt{\langle
R_{g,\perp}^{2}\rangle})$) components, at different degrees of confinement
$\hat{\rm H}$ (see legend in panel (b)). The dashed line in each panel
corresponds to the reference distributions in bulk conditions.

The results for $\langle R_{g}^{2}\rangle$ (Eq. (3)) and the perpendicular and
parallel components, $\langle R_{g,\perp}^{2}\rangle$ and $\langle
R_{g,\parallel}^{2}\rangle$, are reported in Fig. 1. As ${\rm H}$ decreases,
the transverse component $\langle R_{g,\perp}^{2}\rangle$ decreases (green
curve in Fig. 1(a)) as expected. Conversely, the parallel component $\langle
R_{g,\parallel}^{2}\rangle$ grows with confinement (red curve in Fig. 1(a))
because the ring is forced to spread along the plane of the slit. Together,
these two effects produce a characteristic non-monotonic behavior in the
overall $\langle R_{g}^{2}(\hat{\rm H})\rangle$ (blue curve in Fig. 1(a)) with
the minimum attained around $\hat{\rm H}\simeq 0.7$, i.e. where confinement
effects are expected to become more pronounced. Interestingly, for high
confinement ($\hat{\rm H}\lesssim 0.3$), the rings are markedly larger than
the bulk reference (blue dotted curve in Fig. 1(a)). In a previous study
D’Adamo et al. (2017) of randomly concatenated rings under slit confinement,
the non-monotonic behavior was also observed but the swelling compared to the
bulk state was not seen. We attribute this discrepancy to the fact that, in
the previous work, rings without excluded volume were considered which could
have favored more compact conformations.

Beyond average values, we have also computed the corresponding probability
distributions, $P(R_{g})$, $P(R_{g,\perp})$ and $P(R_{g,\parallel})$, and
represented each of them (see Fig. 1, panels (b) to (d)) in the corresponding
scaled variable to ease comparison. While the distributions of the parallel
component of the gyration radius are fundamentally unaffected by confinement
(Fig. 1(c)), the ones of the normal components (see Fig. 1(d)) undergo a
significant change in shape as the confinement becomes stronger, in particular
becoming more peaked. Together these changes produce an interesting effect on
the distributions of the full gyration radius (Fig. 1(b)), which are
characterized by higher tails for the systems under confinement. This suggests
that, under confinement, rings assume more heterogeneous sizes.

Figure 2:  $\langle\Lambda_{1}^{2}\rangle/\langle\Lambda_{1}^{3}\rangle$ and
$\langle\Lambda_{2}^{2}\rangle/\langle\Lambda_{1}^{3}\rangle$, average ring
shapes expressed as the ratios between the two largest eigenvalues of the ring
mean gyration tensor $Q$ (Eq. (2)) and the smallest one, as a function of the
degree of confinement $\hat{\rm H}$ (see Sec. 2.1 for definition). Dotted
lines correspond to the power-law best fits obtained on the first three points
of each curve. Dashed horizontal lines correspond to the bulk reference values
of the two ratios.

We study then rings’ shapes and anisotropies by looking at the ratios: (i)
$\langle\lambda^{2}_{1}\rangle/\langle\lambda^{2}_{2}\rangle$, (ii)
$\langle\lambda^{2}_{1}\rangle/\langle\lambda^{2}_{3}\rangle$ and (iii)
$\langle\lambda^{2}_{2}\rangle/\langle\lambda^{2}_{3}\rangle$. The first ratio
indicates the elongation or “asphericity” of the ring mean shape, while the
other two measure the extent to which rings become effectively flat due to
slit confinement. Results are shown in Fig. 2, where it is clear that at mild
confinement $\hat{\rm H}\gtrsim 1.5$ rings attain the shame shape of the bulk
ones. At higher confinement ($\hat{\rm H}\lesssim 0.7$), the ratios to the
smallest eigenvalues (blue and red curves in Fig. 2) are described by
characteristic power-law behaviors $\sim\\!\hat{\rm H}^{-\alpha}$ with similar
$\alpha$’s: more precisely, the exponent for
$\langle\lambda^{2}_{1}\rangle/\langle\lambda^{2}_{3}\rangle$,
$\alpha=1.831\pm 0.001$, is only slightly larger than the exponent for
$\langle\lambda^{2}_{2}\rangle/\langle\lambda^{2}_{3}\rangle$, $\alpha\simeq
1.781\pm 0.002$. This difference is also evident in the behavior of
$\langle\lambda^{2}_{1}\rangle/\langle\lambda^{2}_{2}\rangle$ (green curve in
Fig. 2), which increases slightly with confinement. In summary, our analysis
shows that rings’ flattening due to confinement biases the chains towards more
elongated shape.

####  III.1.2 3.1.2. Bond-vector correlation function

Figure 3:  (a) $c(\ell)$, bond-vector correlation function as a function of
the contour length distance $\ell$. Colors are for different confinements,
dashed and dot-dashed lines are for bulk melts and melts of non-concatenated
and unknotted rings (see legend). (b) $\langle\cos(\theta)\rangle$, mean
cosine value between two consecutive bonds along the chain as a function of
the degree of confinement $\hat{\rm H}$. (c) $c_{\parallel}(\ell)$,
contribution to the bond-vector correlation function in the $xy$-plane
parallel to the slit. (d) $c_{\perp}(\ell)$, contribution to the bond-vector
correlation function orthogonal to the plane of the slit; in the inset the
same quantity is represented as a function of the ring contour length
normalized by the slit thickness, $\ell/H$. Colors and symbols in panels (c)
and (d) are as in panel (a).

We investigate now in more detail how the folding of polymer chains is
affected by confinement by looking at the bond-vector correlation function,

$c(\ell)\equiv\frac{\langle\vec{t}(\ell^{\prime})\cdot\vec{t}(\ell^{\prime}+\ell)\rangle}{\langle
t(\ell^{\prime})^{2}\rangle}\,,$ (4)

as a function of the polymer contour length $\ell$. This quantity gives useful
insight when applied to bulk $3d$ melts of unknotted and non-concatenated
rings, in particular its distinct Ubertini et al. (2022) anti-correlation is a
symptom of the double folding of the polymer chains at the entanglement scale
(dot-dashed line in Fig. 3(a)). In contrast (dashed line in Fig. 3(a)), bulk
$3d$ melts of randomly knotted and concatenated rings exhibit normal
exponential decay behavior Ubertini and Rosa (2023) and are not characterized
by double folding, hence the anti-correlation is absent.

To investigate the impact of confinement on chain folding, we have computed
$c(\ell)$ for the confined rings. Results (Fig. 3) exhibit several noteworthy
effects. Firstly (Fig. 3(a)), for confined rings at small $\ell$ $c(\ell)$
decays more slowly than the bulk counterpart. This is the consequence (Fig.
3(b)) of the increase of the mean-cosine of the angle between consecutive bond
vectors, $\langle\cos(\theta)\rangle$, as confinements increases: in other
words, confined rings are slightly stiffer than the bulk reference and this
confinement-enhanced stiffness grows with the confinement. At the same time,
$c(\ell)$ develops a characteristic anti-correlation that exhibits non-
monotonic dependence on $\hat{\rm H}$: in particular the deepest minimum
occurs at $\hat{\rm H}\simeq 0.7$, i.e. the same value at which the gyration
radius (Fig. 1(a)) attains its minimum value. Moreover, the minimum itself
disappears at the highest level of confinement. This peculiar behavior can be
explained by considering the individual contributions of the parallel and
transverse components of $c(\ell)$. $c_{\parallel}(\ell)$ does not exhibit any
minima (Fig. 3(c)), while $c_{\perp}(\ell)$ displays a minimum for all values
of $\hat{\rm H}$ (Fig. 3(d)). The mismatch in the values of $\ell$, at which
$c_{\perp}(\ell)$ is minimum while $c_{\parallel}(\ell)\simeq 0$, causes the
non-monotonicity of the full $c(\ell)$. The latter goes to zero for similar
values of $\ell$ for all $\hat{\rm H}$, demonstrating that correlations grow
mildly with the confinement. In contrast, $c_{\perp}(\ell)$ shows a minimum
for $\ell$ close to the thickness of the slit ${\rm H}$ (Fig. 3(d), inset).
This is due to the back-folding of the polymer filaments induced by the
hitting with the impenetrable walls of the slit: of course this effect is more
pronounced under strong confinement conditions, i.e. for ${\rm H}/\ell_{K}\leq
1$. Thus the minima in $c(\ell)$ appear when ${\rm H}$ has similar value to
the correlation length of $c_{\perp}(\ell)$, indicating the competition
between these two length scales.

####  III.1.3 3.1.3. Contact probability

As just shown, confinement alters the metric properties of the polymers. Then,
it is natural to expect that the consequent reorganization of the chains
modifies the intra-chain polymer interactions. To test this hypothesis, we
compute the mean contact probability between two monomers at contour length
separation $\ell=n\langle b\rangle$,

$\langle
p_{c}(\ell)\rangle=\left\langle\frac{1}{N}\sum_{i=1}^{N}\Theta(r_{c}-|{\vec{r}}_{i}-{\vec{r}}_{i+n}|)\right\rangle\,,$
(5)

where $\Theta(x)$ is the Heaviside step function and the “contact distance”
$r_{c}$ is set to the unit lattice size $a$ (notice also that periodic
conditions due to the ring geometry are tacitly assumed in Eq. (5)).

Figure 4:  Mean contact probabilities, $\langle p_{c}\rangle$ (Eq. (5)), as a
function of $\xi=\ell\left(1-\ell/L\right)$ where $\ell$ is the contour length
separation between monomers and $L$ is ring total contour length. Colors are
for different confinements, dashed and dot-dashed lines are for bulk melts and
melts of non-concatenated and unknotted rings (see legend). Inset: local
differential exponent $\equiv\frac{d\log\langle p_{c}\rangle}{d\log\xi}$.

Results are shown in Fig. 4, where $\langle p_{c}\rangle$ is plotted against
the “effective” variable $\xi=\ell(1-\ell/L)$ in order to reduce Rosa and
Everaers (2019) finite-size effect due to the ring geometry. First, one can
notice that in bulk systems, as we let rings perform strand crossings, long-
distance contacts decrease (dashed line) with respect to melts of non-
concatenated and unknotted rings (dot-dashed line). In contrast, confinement
leads to an increase in the tail of the mean contact probability compared to
the bulk reference. Notably, at $\hat{\rm H}=0.30$, the tail’s slope is
slightly less steep than in the non-concatenated state.

To get more insight, it is interesting to look at the exponent controlling the
asymptotic power-law decay, $\langle p_{c}\rangle\simeq\xi^{-\gamma}$ (Fig. 4,
inset). In bulk, strand-crossing rings attain ideal statistics characterized
by $\gamma\simeq 1.5$, as confirmed by our previous findings Ubertini and Rosa
(2021). In contrast, confinement leads to a decrease in $\gamma$ which becomes
close to the same asymptotic value as the non-concatenated state,
$\gamma\simeq 1.15$. Based on mean-field arguments Halverson et al. (2014),
$\gamma=d\nu$ where $d$ is space dimension and $\nu$ is the metric exponent of
the chain relating Doi and Edwards (1986); Rubinstein and Colby (2003) the
chain mean linear size to the number of monomers (i.e., $\langle
R^{2}_{g}\rangle\sim N^{2\nu}$). Strand-crossing rings in bulk exhibit ideal
statistics with $\nu=1/2$ Ubertini and Rosa (2021) and they are characterized
by $\gamma=\frac{3}{2}$ in three dimensions. In confined systems, however, the
rings cannot fold freely in three dimensions, effectively reducing the
dimensionality of the system and resulting in a decrease in $\gamma$.

####  III.1.4 3.1.4. Knots statistics

In our kMC algorithm two filaments from the same chain can cross and this
event may induce the formation of a knot along the chain. Characterization of
knots spectra in confined systems have been addressed so far mostly for
isolated chains Tesi et al. (1994); Micheletti et al. (2011); Micheletti and
Orlandini (2012), while less results are available for confined systems at
melt conditions.

Figure 5:  (a) $P_{\rm knot}(\hat{\rm H})$, ring knotting probability (Eq.
(6)) for knots with $K$ crossings and as a function of the degree of
confinement $\hat{\rm H}$. The dashed line correspond to the value for the
bulk melt. (b) $P_{\rm knot}(\hat{\rm H};K)$, probability of finding a knot
with crossing number $K$. Colors are for different confinements, the dashed is
for bulk melts (see legend). $K=0$ correspond to the unknot and $P_{\rm
knot}(\hat{\rm H};K=0)=1-P_{\rm knot}(\hat{\rm H})$ is its corresponding
probability. Knots with $>12$ crossings cannot be distinguished by Topoly
Dabrowski-Tumanski et al. (2021). Composite knots are knots made up by $2$ or
more irreducible knots.

To fill this gap, we have investigated the occurrence of knots by computing
the Jones polynomial of each ring of our systems and, for simplicity, we
present our results based on the number of irreducible crossings (denoted by
$K$, see Sec. 2.3). Specifically, we have computed the probability, $P_{\rm
knot}(\hat{\rm H};K)$, of finding a knot with $K$ irreducible crossings at
given confinement degree $\hat{\rm H}$ and the cumulative knotting
probability:

$P_{\rm knot}(\hat{\rm H})=\sum_{K=3}^{\infty}P_{\rm knot}(\hat{\rm H};K)\,,$
(6)

which gives the probability that a ring in the melt contains a knot (of any
type). As shown in Fig. 5(a), $P_{\rm knot}(\hat{\rm H})$ grows with the
confinement and reaches the maximum value of $\simeq 0.13$ for the smallest
$\hat{\rm H}$, resulting in an increase of $\simeq 130\%$ compared to bulk
reference (dashed line). Both in bulk and in confinement, the most common knot
type is the simplest one, namely the trefoil knot $3_{1}$. Overall (Fig.
5(b)), more complex knots are much less probable for all $\hat{\rm H}$ values,
yet their abundance increases with confinement, see Fig. 5(b) for $P_{\rm
knot}(\hat{\rm H};K)$ and Fig. S2 in SI for the relative population of knot
types with $K$ crossings. In conclusion, our analysis points out that
confinement enhances the probability of knot formation, yet the overall
occurrence of knots (i.e., $P_{\rm knot}$) remains relatively low ($\lesssim
0.13$).

###  III.2 3.2. Chain-chain correlations

####  III.2.1 3.2.1. Chain neighbours

The increase of the long-range intra-chain contacts seen in Fig. 4 may be
indicative of the fact that confinement reduces the overlap between distinct
chains or, in other words, ring-ring contacts should decrease. To test this
hypothesis, we introduce the variable for the number of neighbors of ring $i$
$(i=1,2,...,M)$,

$\rho_{i}^{\rm ring}\equiv\sum_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{M}\Theta\left(2\sqrt{\langle
R_{g}^{2}\rangle}-|{\vec{r}}_{{\rm CM},i}-{\vec{r}}_{{\rm CM},j}|\right)\,,$
(7)

where $\Theta(x)$ is the Heaviside step function, $\langle R^{2}_{g}\rangle$
is the mean square gyration radius of the system, and ${\vec{r}}_{{\rm CM},j}$
represents the centre of mass position of the $j$-th ring. According to Eq.
(7), two rings are defined as “neighbor” whenever the spatial distance between
their centres of mass is smaller than the twice the root mean-square gyration
radius of the system. We have measured the distribution function of $\rho^{\rm
ring}$, $P(\rho^{\rm ring})$ and its mean value, $\langle\rho^{\rm
ring}\rangle$, at different confinements and we study these quantities in
relation to the distribution of spatial distances between the centres of mass
$d_{\rm CM-CM}$ for neighboring rings, $P(d_{\rm CM-CM}\,|\,{\rm neighbors})$.

Figure 6:  (a) Distribution function, $P(\rho^{\rm ring})$, of the the number
of neighbors per chain $\rho^{\rm ring}$. Inset: mean number of neighbors per
ring, $\langle\rho^{\rm ring}\rangle$. (b) Distribution function of the
distances between the centres of mass of neighbour chains, $P(d_{\rm CM-
CM}\,|\,{\rm neighbours})$, as a function of the variable normalized to twice
the root mean-square gyration radius, $2\sqrt{\langle R_{g}^{2}\rangle}$ (Eq.
(3)), of the rings. Colors are for different confinements, dashed and dot-
dashed lines are for bulk melts and melts of non-concatenated and unknotted
rings (see legend).

Results are shown in Fig. 6, from which it is evident (panel (a)) that
$\langle\rho^{\rm ring}\rangle$ decreases as confinement increases, with
$\langle\rho^{\rm ring}\rangle$ being always smaller with respect to the bulk
reference (dashed line) and even smaller (for the tighter confinements
$\hat{\rm H}\lesssim 1.5$) with respect to the non-concatenated and unknotted
case (dot-dashed line). At the same time (panel (b)), the distributions of
spatial distances $d_{\rm CM-CM}$ demonstrate that neighboring chains tend to
overlap more with each other under stronger confinement. Taken all together,
we can motivate the reason why the inter-chain contacts decrease in terms of
the geometry of the slit. First, confinement can prevent the formation of
stacked conformations along the transverse direction (see Fig. S3 in SI), and
this surely reduces the inter-chain contacts. Moreover, we observe that, by
reducing the width of the slit, inter-ring distances tend to increase and this
is an effect due to the increasing asymmetry of the slit as confinement
increases (see Fig. S4 in SI).

####  III.2.2 3.2.2. Links

The reduction of inter-chain contacts should also have consequences on the
linking properties of the confined systems. To explore this aspect, we adopt
the approach developed by us in Ref. Ubertini and Rosa (2023) and compute,
$\langle n_{\rm 2link}(|{\rm GLN}|)\rangle$, the mean number of two-chain
links at absolute Gauss linking number $|{\rm GLN}|$ and the mean number of
distinct three-chain links, $\langle n_{\rm 3link}\rangle$, with given chain
topology.

Figure 7:  (a) $\langle n_{\rm 2link}(|{\rm GLN}|)\rangle$, mean number of
links per ring with absolute Gauss linking number $|{\rm GLN}|$. (b) $\langle
n_{\rm 3link}\rangle$, mean number of different three-chain linked structures
per ring. Different colors are for the different confinements, the dashed line
is for the bulk system.

Results for $\langle n_{\rm 2link}(|{\rm GLN}|)\rangle$ are summarized in
panel (a) of Fig. 7. We notice that ring-ring links are mostly Hopf-like
(i.e., with $|{\rm GLN}|=1$) and that confinement reduces the extent to which
the rings are linked, in agreement with the decrease of overlaps between
neighboring chains. In general, the participation in more complex links
decreases exponentially but the rate of decay depends on the level of
confinement in the system. Chains under stronger confinement are characterized
by a slower decay, which can be attributed to the fact that neighboring chains
penetrate each other more (see Fig. 6(b)). Additionally, links with $|{\rm
GLN}|=0$ (i.e., the so-called Whitehead links) have been found between those
with $|{\rm GLN}|=2$ and $|{\rm GLN}|=3$ at all confinements. We further
classify these links by computing their Jones polynomial and determining their
relative abundances (panel (a) in Fig. S5 in SI). We found that, even in this
case, rings under stronger confinement form more complex links with greater
ease.

To examine three-chain links, it is necessary to distinguish between two
distinct groups of links: those that can be reduced to two-chain links and
irreducible ones Ubertini and Rosa (2023). The first group include: (a)
poly(3)catenanes, chains made of three rings in which two non-concatenated
rings are connected to a common ring, and (b) triangles, triplets of rings
which are all pairwise concatenated. Both (a) and (b) can be detected via
pairwise linking. Instead, irreducible three-chain links cannot be detected
via pairwise linking and can be further divided in two sub-types: (c)
poly(2)catenane+1-ring, structures made of a poly(2)catenane plus another ring
which is not directly concatenated (in a pairwise manner) to any of the other
two, and (d) Brunnian links, non-trivial links which become a set of trivial
links whenever one component ring is unlinked from the others (the so called
Borromean conformation, the link $6^{3}_{2}$, constitutes the easiest example
of this kind). By resorting to the shrinking method described in Ubertini and
Rosa (2023), we have detected links beloning to the last two classes and
computed $\langle n_{\rm 3link}\rangle$ for the different types of three-chain
links (Fig. 7(b)). It is clear from $\langle n_{\rm 3link}\rangle$ that links
organize onto a network made almost entirely via pairwise concatenation both
in bulk and in confinement. Irreducible three-chain links are much more rare
and decrease with the degree of confinement, for this reason subsequent
analysis has been performed by neglecting these three-chain links
contribution. A detailed topological classification of these structures has
been reported in Fig. S5(b) in SI, and even in this case three-chain links
with higher crossings seem to be more likely for more confined systems.

####  III.2.3 3.2.3. Polymer network and entanglements

Concatenated rings give rise to a fully connected polymer network Bobbili and
Milner (2020); Ubertini and Rosa (2023). To characterize this network, we
define Ubertini and Rosa (2023) the linking degree ${\rm LD}_{i}$ of ring $i$,

${\rm LD}_{i}=\sum_{j=1}^{M}\chi_{ij}\,C_{ij}\,,$ (8)

where the sum runs over the total number of chains in the melt, and where
$C_{ij}$ is the $M\times M$ matrix expressing the concatenation status between
rings $i$ and $j$:

$C_{ij}=\left\\{\begin{array}[]{cl}0\,,&\mbox{if $i=j$}\\\ \\\ 1\,,&\mbox{if
$i\neq j$ {\it and} form a two-chain link}\\\ \\\
0\,,&\mbox{otherwise}\end{array}\right.$ (9)

The “weight” factor $\chi_{ij}$ takes into account the “complexity” of two-
chain links: $\chi_{ij}=|{\rm GLN}|$ or $=\frac{K}{2}$ depending on whether
${\rm GLN}\neq 0$ or ${\rm GLN}=0$ respectively. Here, $K$ is the number of
crossings characterizing the link or, in other words, each crossing of the
link contributes $1/2$ to an entanglement point. This quantity is of special
interest as we have recently showed Ubertini and Rosa (2023) that the mean
value $\langle{\rm LD}\rangle\equiv\langle\frac{1}{M}\sum_{i=1}^{M}{\rm
LD}_{i}\rangle$ is directly connected to the entanglement length of the melt,
$N_{e}$, via the relation $\langle{\rm LD}\rangle=N/N_{e}$ of the systems. To
complement this analysis, we have also computed the distribution of the values
${\rm LD}$ at the single ring level, $P({\rm LD})$, which gives us information
about the heterogeneity of the network.

Figure 8:  (a) Mean linking degree, $\langle{\rm LD}\rangle$, as a function of
the confinement. The horizontal dotted line represents the bulk value. (b)
Distribution functions, $P({\rm LD})$, of the linking degree as a function of
the variable normalized to the corresponding mean value $\langle{\rm
LD}\rangle$. Different colors are for the different confinements, the dashed
line is for the bulk system.

Results are presented in Fig. 8. $\langle{\rm LD}\rangle$ (panel (a))
decreases as a function of the confinement, up to a reduction of $\simeq 60\%$
with respect to bulk conditions. Then, by looking at the distribution
functions (panel (b)) of the linking degree as a function of $X={\rm
LD}/\langle{\rm LD}\rangle$, we see that the curves at mild confinements
display the same behavior for bulk conditions. Conversely, tails become
stronger for more confined systems. This is in agreement with the behaviour
seen for the distribution functions of the sizes of the rings (Fig. 1(b)),
where the tails are higher for stronger confinements. Fluctuations of ring
size may impact on concatenation since smaller rings will be less concatenated
having less possibility to reach other rings, while bigger rings can host more
contacts and consequently more concatenations. To sum up, the resulting
networks of concatenated rings tend to be more heterogeneous as the
confinement become stronger in line with the fluctuations on the rings’ sizes.

##  IV 4\. Discussion and conclusions

Our findings illustrate the impact that slit confinement has on the spatial
structure of randomly concatenated and knotted ring polymers in melt
conditions.

At the single-chain level, our investigation shows that as rings flatten with
increasing confinement they tend to adopt more elongated conformations. At the
same time, rings become slightly more rigid with the confinement, a tendency
captured by the increase of the correlation ($\langle\cos(\theta)\rangle$,
Fig. 3(b)) between consecutive bonds along the chain. We have also
demonstrated that the competition between the Kuhn length of the polymers,
$\ell_{K}$, and the height of the slight, ${\rm H}$, induces a non-monotonous
behavior on the bond-vector correlation function, $c(\ell)$ (see Fig. 3(a, c,
d)). In general, the impact of confinement on ring conformations becomes
particularly pronounced with respect to the formation of long intra-chain
contacts as the slit narrows (see Fig. 4), resulting in more compact rings.
Finally, these changes have significant repercussions on the knotting
probability which increases with the confinement, and for which we register an
increase of $\simeq 130\%$ compared to the bulk value (see Fig. 5(a)).

The effects of slit confinement on the inter-chain statistics is similarly
noteworthy. Specifically, as the level of confinement increases, the average
number of neighbors per ring, $\langle\rho^{\rm ring}\rangle$, experiences a
considerable decrease (see Fig. 6(a)). This is directly connected to the
decrease of the mean linking degree, $\langle{\rm LD}\rangle$, which displays
a total reduction of $\simeq 60\%$ with respect to bulk conditions.
This finding has two interesting implications. First, being $\langle{\rm
LD}\rangle$ directly related to the mean number of entanglement strands per
ring, the decrease of $\langle{\rm LD}\rangle$ means that, at fixed density,
confinement alone may alter the entanglement properties of the system making
$N_{e}$ effectively bigger. This would explains recent findings Kim et al.
(2021a, b) showing that for both, linear chains and rings in two-dimensional
melts, the resulting dynamical quantities display a quite surprising Rouse-
like behavior Doi and Edwards (1986); Rubinstein and Colby (2003) which,
ultimately, points towards the effective irrelevance of entanglement effects
due to inter-chain interactions. It is worth stressing that, although the mean
number of ring-ring concatenations (i.e., $\langle{\rm LD}\rangle$) decreasing
with confinement is not entirely surprisingly (in $2d$ rings can not be
concatenated), the important point to stress here is that the works Kim et al.
(2021a, b) and our analysis here and in Ubertini and Rosa (2023) suggest that
entanglements are indeed well captured Qin and Milner (2014); Milner (2020);
Bobbili and Milner (2020) by two-chain topological links alone. Finally, it is
worth recalling that the elastic plateau modulus $G_{0}$, which quantifies the
stress-strain relationship of polymeric materials, is related to the total
number of entanglement strands of the melt, $G_{0}\propto\frac{NM}{N_{e}}$
Rubinstein and Colby (2003). Then our results imply that, as confinement
grows, the resulting network becomes softer ($G_{0}$ decreases), highlighting
the important role of geometric confinement on the mechanical properties of
the stored polymer network.

##  V Supporting Information

Time mean-square displacement of monomers in the frame of the centre of mass
of the corresponding ring, fractional population of knot types, contour plots
for the joint distribution function of parallel and transverse components of
the distances between the centres of mass of neighboring rings, distribution
functions of the distances between the rings’ centres of mass, fractional
population of two-chain links with ${\rm GLN}=0$.

##  VI Acknowledgement

The authors acknowledge networking support by the COST Action CA17139
(EUTOPIA).

## References

  * Wu et al. (2017) Q. Wu, P. M. Rauscher, X. Lang, R. J. Wojtecki, J. J. De Pablo, M. J. Hore, and S. J. Rowan, Science 358, 1434 (2017).
  * Hart et al. (2021) L. F. Hart, J. E. Hertzog, P. M. Rauscher, B. W. Rawe, M. M. Tranquilli, and S. J. Rowan, Nature Reviews Materials 6, 508 (2021).
  * Krajina et al. (2018) B. A. Krajina, A. Zhu, S. C. Heilshorn, and A. J. Spakowitz, Phys. Rev. Lett. 121, 148001 (2018).
  * Raphaël et al. (1997) E. Raphaël, C. Gay, and P. G. de Gennes, Journal of Statistical Physics 89, 111 (1997).
  * Chen et al. (1995) J. Chen, C. A. Rauch, J. H. White, P. T. Englund, and N. R. Cozzarelli, Cell 80, 61 (1995).
  * Chiarantoni and Micheletti (2022) P. Chiarantoni and C. Micheletti, Macromolecules 55, 4523 (2022).
  * Doi and Edwards (1986) M. Doi and S. F. Edwards, _The Theory of Polymer Dynamics_ (Clarendon, Oxford, 1986).
  * Rubinstein and Colby (2003) M. Rubinstein and R. H. Colby, _Polymer Physics_ (Oxford University Press, New York, 2003).
  * Rauscher et al. (2020a) P. M. Rauscher, K. S. Schweizer, S. J. Rowan, and J. J. de Pablo, The Journal of Chemical Physics 152, 214901 (2020a).
  * Rauscher et al. (2020b) P. M. Rauscher, K. S. Schweizer, S. J. Rowan, and J. J. De Pablo, Macromolecules 53, 3390 (2020b).
  * Lang et al. (2012) M. Lang, J. Fischer, and J.-U. Sommer, Macromolecules 45, 7642 (2012).
  * Lang et al. (2014) M. Lang, J. Fischer, M. Werner, and J.-U. Sommer, Phys. Rev. Lett. 112, 238001 (2014).
  * Ubertini and Rosa (2021) M. A. Ubertini and A. Rosa, Phys. Rev. E 104, 054503 (2021).
  * Ubertini and Rosa (2023) M. A. Ubertini and A. Rosa, Macromolecules 56, 3354 (2023).
  * Soh and Doyle (2021) B. W. Soh and P. S. Doyle, ACS Macro Letters 10, 880 (2021).
  * Ubertini et al. (2022) M. A. Ubertini, J. Smrek, and A. Rosa, Macromolecules 55, 10723 (2022).
  * (17) Notice that the values for $\langle b\rangle$ between the bulk and the confined rings are all very close to each other, with a slight increase as confinement decreases.
  * Kremer and Grest (1990) K. Kremer and G. S. Grest, The Journal of Chemical Physics 92, 5057 (1990).
  * Jones (1985) V. F. R. Jones, Bulletin of the American Mathematical Society 12, 103 (1985).
  * Dabrowski-Tumanski et al. (2021) P. Dabrowski-Tumanski, P. Rubach, W. Niemyska, B. A. Gren, and J. I. Sulkowska, Briefings in Bioinformatics 22, bbaa196 (2021).
  * Rolfsen (2003) D. Rolfsen, _Knots and links_ (AMS Chelsea Publishing, 2003).
  * D’Adamo et al. (2017) G. D’Adamo, E. Orlandini, and C. Micheletti, Macromolecules 50, 1713 (2017).
  * Rosa and Everaers (2019) A. Rosa and R. Everaers, The European Physical Journal E 42, 1 (2019).
  * Halverson et al. (2014) J. D. Halverson, J. Smrek, K. Kremer, and A. Y. Grosberg, Reports on Progress in Physics 77, 022601 (2014).
  * Tesi et al. (1994) M. Tesi, E. J. van Rensburgs, E. Orlandini, and S. Whittington, Journal of Physics A: Mathematical and General 27, 347 (1994).
  * Micheletti et al. (2011) C. Micheletti, D. Marenduzzo, and E. Orlandini, Physics Reports 504, 1 (2011).
  * Micheletti and Orlandini (2012) C. Micheletti and E. Orlandini, Macromolecules 45, 2113 (2012).
  * Bobbili and Milner (2020) S. V. Bobbili and S. T. Milner, Macromolecules 53, 3861 (2020).
  * Kim et al. (2021a) J. Kim, J. M. Kim, and C. Baig, Polymer 213, 123308 (2021a).
  * Kim et al. (2021b) J. Kim, J. M. Kim, and C. Baig, Soft Matter 17, 10703 (2021b).
  * Qin and Milner (2014) J. Qin and S. T. Milner, Macromolecules 47, 6077 (2014).
  * Milner (2020) S. T. Milner, Macromolecules 53, 1314 (2020).
  * Hoste et al. (1998) J. Hoste, M. Thistlethwaite, and J. Weeks, Math. Intelligencer 20, 33 (1998).

– Supporting Information –
Spatial organization of slit-confined melts of ring polymers with non-
conserved topology: A lattice Monte Carlo study

Mattia Alberto Ubertini and Angelo Rosa

Figure S1:  Time mean-square displacement of monomers in the centre of mass of
the corresponding ring Kremer and Grest (1990); Ubertini and Rosa (2021),

$\mbox{\hskip 113.81102pt}g_{2}(\tau_{\rm
MC})\equiv\left\langle\frac{1}{M}\sum_{m=1}^{M}\left[({\vec{r}}_{m}(\tau_{\rm
MC})-{\vec{r}}_{\rm CM}(\tau_{\rm MC}))-({\vec{r}}_{m}(0)-{\vec{r}}_{\rm
CM}(0))\right]^{2}\right\rangle\,,$

as a function of the Monte Carlo (MC) time, $\tau_{\rm MC}$. $\vec{r}_{m}$ is
the spatial position of monomer $m$, while $\vec{r}_{\rm CM}$ is the spatial
position of the centre of mass of the ring. All examined systems display a
plateau, that is indicative of their successful equilibration. The time scale
associated to the onset of the plateau is used to estimate which amount from
the initial portion of the corresponding MC trajectory has to be discarded in
order to compute rings’ properties accurately.  Figure S2:  Fractional
population of each knot type – named according to the Rolfsen convention
Rolfsen (2003) – at fixed number of crossings $K$ and for the different values
of confinement $\hat{\text{H}}$ and in bulk conditions (see legend). For
$K=11$, knots are categorized according the Hoste-Thistlethwaite table Hoste
et al. (1998) that split them in alternating ($K_{a\\_i}$) and non-alternating
ones ($K_{n\\_i}$), with the index $i$ used to enumerate the rings separately
within each group. Large relative error bars are due to the limited size of
the sample (notice that the values of each $y$-axis have to be multiplied by
the power-law reported on the top left corner of the corresponding panel).
Figure S3:  Contour plots for the joint distribution function of the slit-
parallel and slit-transverse (or, slit-perpendicular) components of the
distances between the centres of mass of neighboring rings $P(d_{{\rm CM}-{\rm
CM},\parallel},d_{{\rm CM}-{\rm CM},\perp}\,|\,{\rm neighbors})$. Distances
have been rescaled by the corresponding root-mean-square gyration radius of
the rings.  Figure S4:  Distribution functions of the distances between the
rings’ centres of mass, $P(d_{{\rm CM}-{\rm CM}})$. Figure S5:  (a)
$P(K_{i}^{2}|{\rm GLN}=0)$, fractional population of two-chain links
$K_{i}^{2}$ (termed according to the Rolfsen convention Rolfsen (2003)) having
${\rm GLN}=0$. Labels have been put for links corresponding to peaks of the
distributions. (b) $P(K_{i}^{3}|\mbox{irreducible})$, fractional population of
three-chain links $K_{i}^{3}$ (also termed according to the Rolfsen convention
Rolfsen (2003)) belonging to the poly(2)catenane+1-ring and Brunnian classes
(see main text for details). Topoly fails Dabrowski-Tumanski et al. (2021)
recognizing links with $>9$ crossings, so these lack categorization. In both
panels, vertical dotted lines delimit areas of links at fixed number of
crossings $K$.