# Neuron Growth Control by PDE Backstepping:
Axon Length Regulation by Tubulin Flux Actuation in Soma

Cenk Demir  Shumon Koga  Miroslav Krstic C. Demir and M. Krstic are with the
Department of Mechanical and Aerospace Engineering, U.C. San Diego, 9500
Gilman Drive, La Jolla, CA, 92093-0411<EMAIL_ADDRESS>and<EMAIL_ADDRESS>Koga is with the Department of Electrical and Computer Engineering, UC San
Diego, 9500 Gilman Drive, La Jolla, CA, 92093-0411<EMAIL_ADDRESS>
###### Abstract

In this work, stabilization of an axonal growth in a neuron associated with
the dynamics of tubulin concentration is proposed by designing a boundary
control. The dynamics are given by a parabolic Partial Differential Equation
(PDE) of the tubulin concentration, with a spatial domain of the axon’s length
governed by an Ordinary Differential Equation (ODE) coupled with the tubulin
concentration in the growth cone. We propose a novel backstepping method for
the coupled PDE-ODE dynamics with a moving boundary, and design a control law
for the tubulin concentration flux in the soma. Through employing the Lyapunov
analysis to a nonlinear target system, we prove a local exponential stability
of the closed-loop system under the proposed control law in the spatial
$H_{1}$-norm.

##  I Introduction

Neuroscience has become one of the most significant areas in biology, as
scientists pursue the understanding of the functionality of perception and the
brain, and investigate therapeutics for neurological diseases and injuries
(see [1, 2]). Neurological disorders, such as Parkinson’s disease, spinal cord
injuries, and Alzheimer’s disease, result from a loss of function of neurons
[3], [4], [5]. Particular medical therapeutics, i.e. Chondroitinase ABC
(ChABC), may restore the functionality of neurons (see [6, 7]) by manipulating
the extracellular matrix (ECM), the area surrounding the neuron which contains
various macromolecules and minerals that facilitate cell processes [8].

Neurons are cells that are specialized to obtain perception by transmitting
electrical signals along their axons. These signals enter from the dendrites,
travel through the axon, and transmit through the growth cone to another
neuron as shown in Figure 1. The presence of chemical cues surrounding the
growth cone attract or repel the axon toward another neuron [9]. When the
chemical cues attract the growth cone, tubulin proteins assemble from free
tubulin dimers and create microtubules which extend the axon toward the other
neuron [10]. The formation of microtubules is determined by the size and
dynamics of tubulin concentration in the neuron and is supported by ECM [11,
12]. The dynamics of tubulin depend on the tubulin production rate in the
soma, the degradation rate, the assembly rate, and the transportation process
[9]. According to recent research, it is also possible to control axon
elongation, namely tubulin concentration, by manipulating ECM [13].

The dynamics of axon growth have been described by several different models
using various mathematical tools. One of the pioneer axonal growth models has
been proposed in [14] by considering the transportation of tubulin as
diffusion. A continuum model of axon growth dynamics has been proposed in
[15], and the stability analysis for the proposed model is studied in [16]. In
another work, the tubulin concentration is modelled by a PDE with a moving
boundary [17, 18]. [19] provides a numerical solution for this moving boundary
PDE model.

While PDEs have been utilized for the computational modeling of axon growth,
stabilization of axon growth by means of control theory has not been studied
so far. Boundary control of PDEs has been intensively developed by the method
of “backstepping” in the last few decades for various kinds of systems [20].
One of the initial contributions was achieved in [21] by applying a Volterra
type of transformation to parabolic PDEs by utilizing the method of successive
approximations. Following the procedure, the class of the system has been
extended to a cascade and coupled PDE-ODE system, see [22, 23, 24]. While most
literature on backstepping design of PDEs has focused on a system with a
constant size of domain in time, recently the method has been successfully
applied to the Stefan problem, which is a special class of a parabolic PDE
with a moving boundary, see [25, 26]. For the model related to axon growth
mentioned earlier, the backstepping method for the Stefan problem has been
designed and applied to a screw extrusion process of a polymer 3-D printing
[27]. The results mentioned above have been proven to achieve global stability
by virtue of the monotonicity of the moving boundary.

Several researchers have tackled stability analysis for nonlinear PDE systems
under the backstepping design of a boundary control by restricting the region
of attraction in a local sense. For instance, [28] proved a local exponential
stability of a $2\times 2$ quasi-linear hyperbolic PDE under the backstepping
design by analyzing a Lyapunov function of the spatial $H_{2}$-norm. As a
class of PDEs with a moving boundary, [29] designed a backstepping control law
for $2\times 2$ hyperbolic PDEs with a moving boundary governed by an ODE
modeling a piston movement, and proved a local exponential stability in the
spatial $H_{1}$-norm by applying the Lyapunov method to the target system. A
similar approach has been done for a moving boundary hyperbolic PDE modeling a
shock-wave arising in traffic congestion [30]. However, those results for
local stability analysis have been achieved only for hyperbolic PDE systems,
even though the axon growth model proposed in [18] is a nonlinear parabolic
PDE system.

In this paper, we develop a boundary control for a coupled PDE-ODE system with
a moving boundary which models the dynamics of tubulin concentration and axon
growth. First, we present a steady-state solution of the tubulin concentration
for a given constant axon length, and obtain a reference error system to be
stabilized at zero states. Next, we apply linearization to the reference error
system to deal with the algebraic nonlinearity. Then, a backstepping
transformation is employed to the linearized reference error dynamics. By
solving the gain kernel equations that are derived from backstepping
transformation, the control law is obtained. Finally, we prove local
exponential stability by applying the Lyapunov method to the nonlinear target
system, which ensures the local stability of the original PDE-ODE system of
the axon growth model.

This paper is structured as follows. Section II introduces the PDE-ODE model
of axon growth and tubulin concentration with the steady-state solution.
Section III presents the control design by the method of backstepping, and the
stability result and its proof are given in Section IV. The paper ends with
the conclusion in Section VI.

Figure 1: Neuron Structure and PDE domain

##  II Modeling of Axon Growth

In this section, we present a mathematical model of axon growth governed by a
moving boundary PDE, derive a steady-state solution for a given set point of
the axon length, and provide a reference error system to be stabilized.

###  II-A Axon growth model by a moving boundary PDE

Tubulin is a group of proteins which is responsible for the growth of a newly
created axon. Two assumptions can be described to model this responsibility.
First, tubulin proteins are modeled as a homogeneous continuum because free
tubulin molecules are very small. Then, tubulin molecules are assumed to be
the only factor responsible for axonal growth. With these assumptions, as
proposed in [18, 19], the axonal growth of a newborn axon by tubulin can be
modelled 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),$ (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)

where the tubulin concentration in the axon is denoted as $c(x,t)$, and the
variables with subscripts $t$ and $x$ denote the partial derivatives with
respect to the subscripts. The variable $q$ denotes the concentration flux.
Subscript “s” is used for the soma of the neuron, and subscript “c” is used
for the cone of the neuron. Namely, $q_{\rm s}(t)$ is the concentration flux
of tubulin in the soma, and $c_{\rm c}(t)$ is the concentration of tubulin in
the cone. The length of the axon in $x$-coordinate is $l(t)$. As time passes,
tubulin in the neuron degrade at the constant rate, $g$. $D$ is the
diffusivity constant, and $a$ is the velocity constant of tubulin proteins in
(1). $l_{\rm c}$ is the growth ratio of the cone, and $\tilde{r}_{\rm g}$ is
the reaction rate to create microtubules. $\tilde{s}_{\rm g}$ is the
disassemble rate which means it is the transformation rate from microbules to
tubulin dimers and $r_{\rm g}$ is a lumped parameter. The equilibrium of the
tubulin concentration in the cone, which is denoted as $c_{\infty}$, causes
the axonal growth to stop. The control problem to be solved in this paper is
presented in the following statement.

Problem: Develop a feedback control law of $q_{\rm s}(t)$ so that $l(t)$
converges to $l_{\rm s}$ for a given desired length of the axon $l_{\rm s}>0$,
subject to the governing equations (1)–(5).

###  II-B Steady-state solution

To tackle the problem stated above, we first solve a steady-state solution of
the concentration profile for a given axon length $l_{\rm s}$. By setting the
time derivatives in (1), (II-A), and (5) to zero, one can derive the steady-
state solution of (1)-(5) as

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

where

$\displaystyle\lambda_{+}=$
$\displaystyle\frac{a+\sqrt{a^{2}+4Dg}}{2D},\quad\lambda_{-}=\frac{a-\sqrt{a^{2}+4Dg}}{2D},$
(7) $\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}}.$ (8)

The steady-state input for the concentration flux in the soma is obtained as

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

###  II-C Reference error system

Let $u(x,t)$, $z_{1}(t)$, and $z_{2}(t)$ be the reference error states, and
$U(t)$ be the reference error input, defined as

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

By substituting the steady-state solution in (6) from the governing equations
(1)-(5), the reference error system is obtained as in [18, 26].

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}.$ (14)

Applying the linearization of $X(t)$ around zero states leads to the following
linearized reference error system (see Section 12-2 in [26] for the detailed
derivation):

$\displaystyle u_{t}(x,t)=$ $\displaystyle Du_{xx}(x,t)-au_{x}(x,t)-gu(x,t),$
(15) $\displaystyle u_{x}(0,t)=$ $\displaystyle U(t),$ (16) $\displaystyle
u(l(t),t)=$ $\displaystyle C^{\top}X(t),$ (17) $\displaystyle\dot{X}(t)=$
$\displaystyle AX(t)+Bu_{x}(l(t),t),$ (18)

where

$\displaystyle A=$ $\displaystyle\left[\begin{array}[]{cc}\tilde{a}&0\\\
r_{\rm g}&0\end{array}\right],\quad B=\left[\begin{array}[]{c}-\beta\\\
0\end{array}\right],$ (23) $\displaystyle C=$
$\displaystyle\left[1\quad-\frac{(a-gl_{\rm c})c_{\infty}}{D}\right]^{\top}.$
(24)

##  III Control Design

The control design in this paper is based on a backstepping transformation
[20], which maps the reference error system $(u,X)$ to a target system
$(w,X)$. The backstepping transformation and the target system in this paper
are given in the remainder of this section.

###  III-A Backstepping transformation and target system

Following the procedure in [25] for the Stefan problem, we consider the
following backstepping transformation:

$\displaystyle w(x,t)=$ $\displaystyle
u(x,t)-\int_{x}^{l(t)}k(x,y)u(y,t)dy-\phi(x-l(t))^{\top}X(t),$ (25)

where $k(x,y)\in{\mathbb{R}}$ and $\phi(x-l(t))\in{\mathbb{R}}^{2}$ are the
gain kernel functions to be determined. We suppose the desired target system
as

$\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)),$ (26) $\displaystyle w_{x}(0,t)=$
$\displaystyle\gamma w(0,t),$ (27) $\displaystyle w(l(t),t)=$ $\displaystyle
0,$ (28) $\displaystyle\dot{X}(t)=$
$\displaystyle(A+BK^{\top})X(t)+Bw_{x}(l(t),t),$ (29)

where $K\in{\mathbb{R}}^{2}$ is a feedback control gain vector chosen to make
$A+BK$ Hurwitz which means that $(A,B)$ is controllable. In detail, by setting

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

one can observe that the conditions for $(k_{1},k_{2})$ makes $A+BK$ Hurwitz.
Also, the redundant nonlinear term $F(x,X(t))\in{\mathbb{R}}$ in (III-A)
caused by the moving boundary is described as

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

###  III-B Gain kernel solution

Taking the time and spatial derivatives of (25) together with the solution of
(15)-(18), and substituting $x=l(t)$ in both the transformation (25) and its
spatial derivative, and by matching with the target system (III-A)–(29), we
have the following PDE and an ODE for gain kernels.

$\displaystyle
k_{xx}(x,y)-k_{yy}(x,y)=\frac{a}{D}\left(k_{x}(x,y)+k_{y}(x,y)\right),$ (32)
$\displaystyle k_{x}(x,x)+k_{y}(x,x)=0,$ (33) $\displaystyle
k(x,l(t))=-\frac{1}{D}\phi(x-l(t))^{\top}B,$ (34) $\displaystyle
D\phi^{\prime\prime}(x-l(t))^{\top}-a\phi^{\prime}(x-l(t))^{\top}-\phi(x-l(t))^{T}\left[gI+A\right]$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \
-Dk_{y}(x,l(t))C^{\top}+ak(x,l(t))C^{\top}=0,$ (35)
$\displaystyle\phi(0)=C,\quad\phi^{\prime}(0)=k(l(t),l(t))C^{\top}+K^{\top}.$
(36)

By the conditions (32)–(34), the solution of $k(x,y)$ is uniquely given by

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

Substituting (37) into (35)–(36), the ODE of $\phi(\cdot)$ becomes

$\displaystyle
D\phi^{\prime\prime}(x-l(t))^{\top}-\phi^{\prime}(x-l(t))^{\top}\left(BC^{\top}+aI\right)$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
-\phi(x-l(t))^{\top}\left[gI+A+\frac{a}{D}BC^{\top}\right]=0,$ (38)
$\displaystyle\phi(0)=C,\quad\phi^{\prime}(0)^{\top}=-\frac{1}{D}C^{\top}BC^{\top}+K^{\top}.$
(39)

The solution to (38)-(39) is given by (see [24])

$\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},$ (40)

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}.$ (41)

In order to check invertibility of (25), we define the inverse Volterra
Transformation

$\displaystyle u(x,t)=$ $\displaystyle
w(x,t)+\int_{x}^{l(t)}q(x,y)w(y,t)dy+\varphi(x-l(t))^{\top}X(t).$ (42)

Then, we apply the same strategy to obtain the inverse gain kernel functions,
$q(x,y)\in\mathbb{R}$, and $\varphi(x-l(t))\in\mathbb{R}^{2}$. The PDE and ODE
for the inverse gain kernels are

$\displaystyle
q_{xx}(x,y)-q_{yy}(x,y)=\frac{a}{D}\left(q_{x}(x,y)+q_{y}(x,y)\right)$ (43)
$\displaystyle q_{x}(x,x)+q_{y}(y,y)=0$ (44) $\displaystyle
q(x,l(t))=-\frac{1}{D}\varphi(x-l(t))^{\top}B$ (45) $\displaystyle
D\varphi^{{}^{\prime\prime}}(x-l(t))^{\top}+a\varphi^{{}^{\prime}}(x-l(t))^{\top}$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \
+\left(gI+A+BK^{\top}\right)\varphi(x-l(t))^{\top}=0$ (46)
$\displaystyle\varphi(0)=C,\quad\varphi^{\prime}(0)=K$ (47)

which is well-posed, so one can obtain the solution of (43)-(47) by applying
the same procedure as the one we applied to the forward kernel equations.

###  III-C Control law

The control design is derived from the boundary condition (27) of the target
system at $x=0$. Substituting $x=0$ into the transformation (25) and its
spatial derivative, and using boundary conditions for both the target system
and the error system, the control input is described as follows

$\displaystyle U(t)=$
$\displaystyle\left(\gamma-\frac{\beta}{D}\right)u(0,t)-\frac{1}{D}\int_{0}^{l(t)}p(x)Bu(x,t)dx$
$\displaystyle+p(l(t))X(t),$ (48)

where $p(x)=\phi^{\prime}(-x)^{\top}-\gamma\phi(-x)^{\top}$. Plugging the
system matrices $(A,B,C,K)$ into the gain kernel function (40), and
calculating the transition matrix, one can explicitly derive the function
$p(x)\in{\mathbb{R}}^{2}$.

##  IV Proof of Main Result

This section presents the main result of the paper, which is a local stability
of the closed-loop system, with providing its proof by considering the
$H_{1}$-norm

$\displaystyle Z(t)=||u(\cdot,t)||^{2}+||u_{x}(\cdot,t)||^{2}+X^{\top}X.$ (49)

We present our main theorem below.

######  Theorem 1.

Consider the closed-loop system consisting of the plant (15)–(18) with the
control law (48) and (40). Suppose the control parameter $\gamma>0$ is chosen
to satisfy $\gamma\geq\frac{a}{D}$. Then, there exist positive parameters
$\bar{M}>0$, $c>0$, and $\kappa>0$, such that if $Z(0)<\bar{M}$ then the
following norm estimate holds

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

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

###  IV-A Basic idea of proof of Theorem 1

Since the backstepping transformation (25) is invertible, the stability
property of the target system (III-A)-(29) is equivalent to the stability
property of the closed-loop system consisting of the plant (15)-(18) with the
control law (48). The local stability of the target system is studied under
the assumption that the following two conditions hold:

$\displaystyle 0<l(t)\leq\bar{l},$ (51)
$\displaystyle|\dot{l}(t)|\leq\bar{v},$ (52)

for some $\bar{l}>l_{\rm s}>0$ and $\bar{v}>0$. The restricted initial state
will be given later in order to satisfy these conditions for all $t\geq 0$.

###  IV-B Useful inequalities

The following inequalities are used in Lyapunov analysis. Under the condition
(51), and by boundary conditions (27) and (28), Poincare’s inequality is
provided as

$\displaystyle||w||^{2}\leq$ $\displaystyle
4\bar{l}^{2}||w_{x}||^{2},\quad||w_{x}||^{2}\leq$ $\displaystyle
2\bar{l}\gamma^{2}w(0,t)^{2}+4\bar{l}^{2}||w_{xx}||^{2},$ (53)

and the Agmon’s inequality is given as

$\displaystyle w_{x}(l(t),t)^{2}\leq$ $\displaystyle
2\gamma^{2}w(0,t)^{2}+4\bar{l}||w_{xx}||^{2}.$ (54)

###  IV-C Proof of Lyapunov stability

We consider the Lyapunov function of the target system as

$\displaystyle V=d_{1}V_{1}+V_{2}+\frac{\gamma}{2}w(0,t)^{2}+d_{2}V_{3},$ (55)

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

$\displaystyle V_{1}=$
$\displaystyle\frac{1}{2}||w||^{2}:=\frac{1}{2}\int_{0}^{l(t)}w(x,t)^{2}dx,$
(56) $\displaystyle V_{2}=$
$\displaystyle\frac{1}{2}||w_{x}||^{2}:=\frac{1}{2}\int_{0}^{l(t)}w_{x}(x,t)^{2}dx,$
(57) $\displaystyle V_{3}=$ $\displaystyle X(t)^{\top}PX(t),$ (58)

where positive definite matrix $P$ is the solution of the following Lyapunov
equation

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

for some positive definite matrix $Q$. Since $(A+BK^{\top})$ is Hurwitz,
positive definite matrices $P$ and $Q$ exist. Due to the positive
definiteness, it holds that

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

where $\lambda_{\rm min}(P)>0$ and $\lambda_{\rm max}(P)>0$ are the smallest
and the largest eigenvalues of $P>0$. Then we have the following lemma.

######  Lemma 1.

Assume that (51)–(52) are satisfied with

$\displaystyle\bar{v}=\min\left\\{\frac{g}{4\gamma},\frac{D}{8\bar{l}}\right\\},$
(61)

for all $t\geq 0$. Then, for sufficiently large $d_{1}>0$ and small $d_{2}>0$,
there exists a positive constant $\beta>0$ such the following norm estimate
holds for all $t\geq 0$:

$\displaystyle\dot{V}\leq-\alpha V+\beta V^{3/2},$ (62)

where
$\alpha=\min\left\\{2g+\frac{D}{4\bar{l}},\frac{4g+d_{1}D}{2},\frac{\lambda_{\rm
min}(Q)}{2\lambda_{\rm
max}(P)},\frac{d_{2}\left(2d_{1}D+g\right)}{4}\right\\}$.

###### Proof.

Taking the time derivative of the Lyapunov functions along the target system
(III-A)–(29), we have

$\displaystyle\dot{V}_{1}=$
$\displaystyle-D||w_{x}||^{2}-g||w||^{2}-\left(\gamma
D-\frac{a}{2}\right)w(0,t)^{2}$
$\displaystyle+\dot{l}(t)\int_{0}^{l(t)}F(x,X(t))w(x,t)dx.$ (63)
$\displaystyle\dot{V}_{2}=$
$\displaystyle-D||w_{xx}||^{2}+a\int_{0}^{l(t)}w_{xx}(x,t)w_{x}(x,t)dx$
$\displaystyle-g||w_{x}||^{2}-\gamma
w(0,t)w_{t}(0,t)-\frac{1}{2}\dot{l}(t)w_{x}(l(t),t)^{2}$
$\displaystyle-\dot{l}(t)(F(l(t),X(t))w_{x}(l(t),t)-\gamma F(0,X(t))w(0,t))$
$\displaystyle-$ $\displaystyle\gamma
gw(0,t)^{2}-\dot{l}(t)\int_{0}^{l(t)}F_{x}(x,X(t))w_{x}(x,t)dx$ (64)
$\displaystyle\dot{V}_{3}=$
$\displaystyle-X(t)^{\top}QX(t)+2w_{x}(l(t),t)B^{\top}PX(t).$ (65)

Applying the Agmon’s inequality (54) and Young’s inequality to $\dot{V}_{2}$
in (64) and $\dot{V}_{3}$ in (65) leads to

$\displaystyle\dot{V}_{2}\leq$
$\displaystyle-\frac{D}{4}||w_{xx}||^{2}-\frac{\gamma
g}{2}w(0,t)^{2}-\left(g-\frac{a^{2}}{D}\right)||w_{x}||^{2}$
$\displaystyle+\dot{l}(t)\gamma F(0,X(t))w(0,t)$
$\displaystyle+\dot{l}(t)\int_{0}^{l(t)}F_{x}(x,X(t))w_{x}(x,t)dx$
$\displaystyle+\frac{\big{|}\dot{l}(t)\big{|}}{2}F(l(t),X(t))^{2}-\gamma
w(0,t)w_{t}(0,t),$ (66) $\displaystyle\dot{V}_{3}\leq$
$\displaystyle-\frac{\lambda_{\rm min}(Q)}{2}X^{\top}X$
$\displaystyle+\frac{2\big{|}B^{\top}P\big{|}^{2}}{\lambda_{\rm
min}(Q)}\left(2\gamma^{2}w(0,t)^{2}+4\bar{l}||w_{xx}||^{2}\right).$ (67)

By using (63), (66) and (67), recalling $\gamma\geq\frac{a}{D}$, and choosing
the constants $d_{1}$ and $d_{2}$ to satisfy

$\displaystyle d_{1}\geq\frac{2a^{2}}{D^{2}},\
d_{2}\leq\min\left\\{\frac{D\lambda_{\rm
min}(Q)}{64\bar{l}\big{|}B^{T}P\big{|}^{2}},\frac{gD\lambda_{\rm
min}(Q)}{16a\left|B^{T}P\right|^{2}}\right\\},$ (68)

the time derivative of Lyapunov function (55) for the target system is shown
to satisfy the following inequality

$\displaystyle\dot{V}$ $\displaystyle\leq-\alpha
V+\frac{\big{|}\dot{l}(t)\big{|}}{2}F(l(t),X(t))^{2}+\dot{l}(t)\gamma
F(0,X(t))w(0,t)$
$\displaystyle+\dot{l}(t)\int_{0}^{l(t)}F_{x}(x,X(t))w_{x}(x,t)+d_{1}F(x,X(t))w(x,t)dx.$
(69)

Taking the square of (31), it follows that the redundant nonlinear terms in
(69) can be bounded by a quadratic norm of the ODE state. Namely, there exist
positive constants $L_{1}>0$, $L_{2}>0$, and $L_{3}$ such that
$F(0,X(t))^{2}\leq L_{1}X^{\top}X$, $F(l(t),X(t))^{2}\leq L_{1}X^{\top}X$,
$\int_{0}^{l(t)}F_{x}(x,X(t))^{2}dx\leq L_{2}X^{\top}X$,
$\int_{0}^{l(t)}F(x,X(t))^{2}dx\leq L_{3}X^{\top}X$. In addition, $\dot{l}(t)$
can be rewritten as $\dot{l}(t)=r_{\rm g}e_{1}X$. By using these relations and
applying Cauchy-Schwarz inequality and Young’s inequality, one can show that

$\displaystyle\dot{V}$ $\displaystyle\leq-\alpha V+r_{\rm
g}\left(\frac{L_{1}(1+\gamma)+1}{2d_{2}\lambda_{\rm
min}(P)}+(1+L_{3}+d_{1}L_{2})\right)V^{3/2}.$ (70)

Thus, there exists $\beta>0$ such that Lemma 1 holds. ∎

###  IV-D Guaranteeing the conditions for all time

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

######  Lemma 2.

There exists a positive constant $M_{1}>0$ such that in the region
$\Omega_{1}:=\\{(w,X)\in H_{1}\times{\mathbb{R}}^{2}|V(t)<M_{1}\\}$ the
conditions (51) and (52) are satisfied.

###### Proof.

By using (14) and (5), $X(t)$ can be described as
$X(t)=\left[\begin{array}[]{cc}\frac{\dot{l}(t)}{r_{\rm g}}&l(t)-l_{\rm
s}\end{array}\right]^{\top}$. 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.$ (71)

The former inequality tells that if $r<\frac{\bar{v}}{r_{\rm g}}$ then the
condition (51) holds. Moreover, the latter 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 the condition (52) holds. Therefore, the constant $r>0$ is chosen
as $r=\min\left\\{\frac{\bar{v}}{r_{\rm g}},l_{\rm s},\bar{l}-l_{\rm
s}\right\\}.$ Since $|X|^{2}\leq\frac{1}{\lambda_{\rm
min}(P)}X^{\top}PX\leq\frac{d_{2}}{\lambda_{\rm min}(P)}V$, we deduce that by
setting $M_{1}=\frac{\lambda_{\rm min}(P)}{d_{2}}r^{2}$, if $V(t)<M_{1}$ holds
then $|X|<r$ and thus the conditions (51) and (52) are satisfied, by which we
can conclude Lemma 2. ∎

######  Lemma 3.

There exists a positive constant $M>0$ such that if $V(0)<M$ then the
conditions (51) and (52) are satisfied and the following norm estimate holds:

$\displaystyle V(t)\leq V(0)\exp\left(-\frac{\alpha}{2}t\right).$ (72)

###### Proof.

For a positive constant $M>0$, let $\Omega:=\\{(w,X)\in
H_{1}\times{\mathbb{R}}^{2}|V(t)<M\\}$. By Lemma 2, it is easily shown that if
$M\leq M_{1}$ then $\Omega\subset\Omega_{1}$, and thus the conditions (51) and
(52) are satisfied in the region $\Omega$. Thus, by Lemma 1, the norm estimate
(62) holds. Moreover, by setting $M\leq\frac{\alpha^{2}}{4\beta^{2}}$, we can
see that applying $V(t)<M$ to (62) leads to

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

by which the norm estimate (72) is deduced. Since (72) is a monotonically
decreasing function in time, by setting
$M=\min\\{M_{1},\frac{\alpha^{2}}{4\beta^{2}}\\}$, the region $\Omega$ is
shown to be an invariant set. Thus, if $V(0)<M$, then $V(t)<M$ for all $t\geq
0$, and one can conclude with Lemma 3. ∎

Finally, owing to Lemma 3, and the equivalent norm estimate in the
$H_{1}$-norm between the target and the closed-loop system, the local
stability of the closed-loop system is proved, which completes the proof of
Theorem 1.

##  V Numerical Simulation

(a)

(b)

Figure 2: Closed-loop response under the designed control law.  TABLE I: Biological constants and control parameters Parameter | Value | Parameter | Value
---|---|---|---
$D$ | $10\times 10^{-6}m^{2}/s$ | $\gamma$ | $10^{4}$
$a$ | $1\times 10^{-8}m/s$ | $l_{\rm c}$ | $4\mu m$
$g$ | $5\times 10^{-7}\ s^{-1}$ | $l_{s}$ | $12\mu m$
$r_{\rm g}$ | $1.783\times 10^{-5}\ m^{4}/(mols)$ | $l_{0}$ | $1\mu m$
$c_{\infty}$ | $0.0119\ mol/m^{3}$ | $k_{1}$ | $-0.1$
$\tilde{r}_{\rm g}$ | $0.053$ | $k_{2}$ | $10^{13}$

Simulation is performed for the axon growth dynamics (1)-(5) under the
designed control law (48). We use the biological constants given in [18], and
choose the control parameters, as shown in Table I. The initial tubulin
concentration is set as a constant profile $c_{0}(x)=2c_{\infty}$, and the
initial axon length is set as $l_{0}=1\mu m$.

Fig. LABEL:subfig-2 shows that the axon length converges to the desired length
$l_{\rm s}$. In addition, Fig. LABEL:subfig-1 depicts that the tubulin
concentration converges to the steady-state solution. Hence, we observe that
the simulation result is consistent with our theoretical results.

##  VI Conclusion

In this paper, a boundary feedback control for an axonal growth model governed
by a coupled PDE-ODE with a moving boundary is studied. The backstepping
transformation is utilized for the original system to convert it to the target
system which has a stable system matrix in ODE. The gain kernel solutions in
the transformation are obtained, and the boundary feedback control law is
designed explicitly. Using Lyapunov’s method, we first prove that the target
system is locally exponentially stable in the spatial $H_{1}$-norm and then we
prove local stability of the original coupled PDE-ODE system of axonal growth
model. Showing the local stability without applying linearization will be
studied in future.

## References

  * [1] L. Ribar and R. Sepulchre, “Neuromorphic control,” _arXiv preprint arXiv:2011.04441_ , 2020.
  * [2] G. Drion, T. O’Leary, J. Dethier, A. Franci, and R. Sepulchre, “Neuronal behaviors: A control perspective,” in _2015 54th IEEE Conference on Decision and Control (CDC)_ , 2015, pp. 1923–1944.
  * [3] R. B. Maccioni, J. P. Muñoz, and L. Barbeito, “The molecular bases of alzheimer’s disease and other neurodegenerative disorders,” _Archives of medical research_ , vol. 32, no. 5, pp. 367–381, 2001.
  * [4] W. Dauer and S. Przedborski, “Parkinson’s disease: mechanisms and models,” _Neuron_ , vol. 39, no. 6, pp. 889–909, 2003.
  * [5] X. Z. Liu, X. M. Xu, R. Hu, C. Du, S. X. Zhang, J. W. McDonald, H. X. Dong, Y. J. Wu, G. S. Fan, M. F. Jacquin, _et al._ , “Neuronal and glial apoptosis after traumatic spinal cord injury,” _Journal of Neuroscience_ , vol. 17, no. 14, pp. 5395–5406, 1997.
  * [6] S. Karimi-Abdolrezaee, E. Eftekharpour, J. Wang, D. Schut, and M. G. Fehlings, “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_ , vol. 30, no. 5, pp. 1657–1676, 2010.
  * [7] E. J. Bradbury and L. M. Carter, “Manipulating the glial scar: chondroitinase abc as a therapy for spinal cord injury,” _Brain research bulletin_ , vol. 84, no. 4-5, pp. 306–316, 2011.
  * [8] C. Frantz, K. M. Stewart, and V. M. Weaver, “The extracellular matrix at a glance,” _Journal of cell science_ , vol. 123, no. 24, pp. 4195–4200, 2010\.
  * [9] S. Diehl, E. Henningsson, A. Heyden, and S. Perna, “A one-dimensional moving-boundary model for tubulin-driven axonal growth,” _Journal of Theoretical Biology_ , vol. 358, pp. 194–207, 2014.
  * [10] J.-P. Julien, “Neurofilament functions in health and disease,” _Current opinion in neurobiology_ , vol. 9, no. 5, pp. 554–560, 1999.
  * [11] E. W. Dent and F. B. Gertler, “Cytoskeletal dynamics and transport in growth cone motility and axon guidance,” _Neuron_ , vol. 40, no. 2, pp. 209–227, 2003.
  * [12] C. S. Barros, S. J. Franco, and U. Müller, “Extracellular matrix: functions in the nervous system,” _Cold Spring Harbor perspectives in biology_ , vol. 3, no. 1, p. a005108, 2011.
  * [13] E. Burnside and E. Bradbury, “Manipulating the extracellular matrix and its role in brain and spinal cord plasticity and repair,” _Neuropathology and applied neurobiology_ , vol. 40, no. 1, pp. 26–59, 2014.
  * [14] M. P. Van Veen and J. Van Pelt, “Neuritic growth rate described by modeling microtubule dynamics,” _Bulletin of mathematical biology_ , vol. 56, no. 2, pp. 249–273, 1994.
  * [15] D. R. McLean, A. van Ooyen, and B. P. Graham, “Continuum model for tubulin-driven neurite elongation,” _Neurocomputing_ , vol. 58, pp. 511–516, 2004.
  * [16] D. R. McLean and B. P. Graham, “Stability in a mathematical model of neurite elongation,” _Mathematical medicine and biology: a journal of the IMA_ , vol. 23, no. 2, pp. 101–117, 2006.
  * [17] B. P. Graham and A. Van Ooyen, “Mathematical modelling and numerical simulation of the morphological development of neurons,” _BMC neuroscience_ , vol. 7, no. 1, pp. 1–12, 2006.
  * [18] S. Diehl, E. Henningsson, A. Heyden, and S. Perna, “A one-dimensional moving-boundary model for tubulin-driven axonal growth,” _Journal of theoretical biology_ , vol. 358, pp. 194–207, 2014.
  * [19] S. Diehl, E. Henningsson, and A. Heyden, “Efficient simulations of tubulin-driven axonal growth,” _Journal of computational neuroscience_ , vol. 41, no. 1, pp. 45–63, 2016.
  * [20] M. Krstic and A. Smyshlyaev, _Boundary control of PDEs: A course on backstepping designs_. SIAM, 2008.
  * [21] A. Smyshlyaev and M. Krstic, “Closed-form boundary state feedbacks for a class of 1-d partial integro-differential equations,” _IEEE Transactions on Automatic Control_ , vol. 49, no. 12, pp. 2185–2202, 2004.
  * [22] M. Krstic, “Compensating actuator and sensor dynamics governed by diffusion pdes,” _Systems & Control Letters_, vol. 58, no. 5, pp. 372–377, 2009\.
  * [23] G. A. Susto and M. Krstic, “Control of pde–ode cascades with neumann interconnections,” _Journal of the Franklin Institute_ , vol. 347, no. 1, pp. 284–314, 2010.
  * [24] S. Tang and C. Xie, “State and output feedback boundary control for a coupled pde–ode system,” _Systems & Control Letters_, vol. 60, no. 8, pp. 540–545, 2011.
  * [25] S. Koga, M. Diagne, and M. Krstic, “Control and state estimation of the one-phase stefan problem via backstepping design,” _IEEE Transactions on Automatic Control_ , vol. 64, no. 2, pp. 510–525, 2018.
  * [26] S. Koga and M. Krstic, _Materials Phase Change PDE Control and Estimation: From Additive Manufacturing to Polar Ice_. Springer Nature, 2020.
  * [27] S. Koga, D. Straub, M. Diagne, and M. Krstic, “Stabilization of filament production rate for screw extrusion-based polymer three-dimensional-printing,” _Journal of Dynamic Systems, Measurement, and Control_ , vol. 142, no. 3, 2020.
  * [28] J.-M. Coron, R. Vazquez, M. Krstic, and G. Bastin, “Local exponential h^2 stabilization of a 2$\backslash$times2 quasilinear hyperbolic system using backstepping,” _SIAM Journal on Control and Optimization_ , vol. 51, no. 3, pp. 2005–2035, 2013.
  * [29] M. Buisson-Fenet, S. Koga, and M. Krstic, “Control of piston position in inviscid gas by bilateral boundary actuation,” in _2018 IEEE Conference on Decision and Control (CDC)_. IEEE, 2018, pp. 5622–5627.
  * [30] H. Yu, M. Diagne, L. Zhang, and M. Krstic, “Bilateral boundary control of moving shockwave in lwr model of congested traffic,” _IEEE Transactions on Automatic Control_ , 2020.