# Dynamic network congestion pricing based on
deep reinforcement learning*

Kimihiro Sato1, Toru Seo2 and Takashi Fuse1 *This work was mainly supported by
a JSPS KAKENHI Grant-in-Aid for Scientific Research 20H02267.1K. Sato and T.
Fuse are with Department of Civil Engineering, the University of Tokyo, Tokyo
113-8656, Japan.2T. Seo is with the Department of Civil and Environmental
Engineering, Tokyo Institute of Technology, Tokyo 152-8552, Japan. Email:
<EMAIL_ADDRESS>
###### Abstract

Traffic congestion is a serious problem in urban areas. Dynamic congestion
pricing is one of the useful schemes to eliminate traffic congestion in
strategic scale. However, in the reality, an optimal dynamic congestion
pricing is very difficult or impossible to determine theoretically, because
road networks are usually large and complicated, and behavior of road users is
uncertain. To account for this challenge, this work proposes a dynamic
congestion pricing method using deep reinforcement learning (DRL). It is
designed to eliminate traffic congestion based on observable data in general
large-scale road networks, by leveraging the data-driven nature of deep
reinforcement learning. One of the novel elements of the proposed method is
the distributed and cooperative learning scheme. Specifically, the DRL is
implemented by a spatial-temporally distributed manner, and cooperation among
DRL agents is established by novel techniques we call spatially shared reward
and temporally switching learning. It enables fast and computationally
efficient learning in large-scale networks. The numerical experiments using
Sioux Falls Network showed that the proposed method works well thanks to the
novel learning scheme.

##  I INTRODUCTION

Traffic congestion is a serious problem in urban areas. It damages the economy
due to the wasted travel time, as well as traffic safety and environment.
Appropriate control schemes that reduce or eliminate congestion are required.

Dynamic congestion pricing [1] is one of the useful schemes to eliminate
traffic congestion in strategic scale. Congestion is mainly caused by spatial
and/or temporal concentration of traffic demand. In a dynamic congestion
pricing scheme, appropriate tolls are charged to road users in certain roads
in certain time duration to mitigate the demand concentration. In the
literature, it has been theoretically shown that an optimal dynamic congestion
pricing could eliminate traffic congestion, if the road network shape is
simple and the behavioral principle of road users is known [1].

In the reality, an optimal dynamic congestion pricing is very difficult or
impossible to determine theoretically, because road networks are usually large
and complicated, and behavior of road users is uncertain. To account for this
challenge, trial-and-error pricing schemes have been proposed [2] [3] [4].
They charge a trial toll first, and then adjust a toll in the next term based
on the road user’s response to the trial toll. By repeating this trial-and-
error process for a certain duration, they eventually find the optimal toll.
Unfortunately, this process sometimes takes a long time (e.g., more than a
year) in large-scale networks, which means that conventional trial-and-error
method may be unpractical in such situations.

The deep reinforcement learning (DRL) can be considered as an efficient data-
driven approach to execute such trial-and-error process [5]. The notable
advantages of DRL are as follows. First, it can automatically extract
meaningful patterns from massive observation data. Second, it can
automatically optimize the adjustment process using the data. These advantages
are useful for dynamic congestion pricing optimization in a large road
network.

The goal of this work is to develop a dynamic congestion pricing method using
DRL. It determines efficient congestion tolls based only on observation data
(e.g., travel time) in general large-scale networks in a day-to-day setting.
We consider a general transportation model in which travelers can choose their
departure time and routes simultaneously. To achieve the goal, this work
develops a novel distributed and cooperative learning scheme for DRL is
developed. Specifically, the DRL is implemented by a spatial-temporally
distributed manner, and cooperation among DRL agents is established by novel
techniques we call spatially shared reward and temporally switching learning.
Thanks to this scheme, the method’s learning speed is computationally
efficient and fast, even in large-scale applications. The proposed method is
validated using actual road network data and compared with other methods, and
the effectiveness of the proposed learning scheme was confirmed.

##  II LITERATURE REVIEW

Application of reinforcement learning (RL) to transportation problems can be
reviewed as follows. Wang et al. [6] apply deep reinforcement learning to
traffic signal control. Their proposed method uses ”high-resolution event-
based data” , which keep track of passage and presence of vehicles by
recording activation/deactivation events of vehicle detectors. Li et al. [7]
proposed ”Knowledge Sharing Deep Deterministic Policy Gradient (KS-DDPG)”,
which introduces a communication protocol for knowledge sharing among multi-
agents, and applies it to cooperative control of traffic signals. Li et al.
[8] propose a DDPG-based driving strategy for individual vehicles to mitigate
oscillations in stop-and-go-waves and optimize traffic safety. It is shown to
reduce collision risk.

Some studies use RL or neural networks to determine dynamic congestion tolls.
Zhu and Ukkusuri [9] apply RL to dynamic congestion pricing with distance-
based tolls. Toll is got closer to the optimal one by progressing a learning
process with traffic flow data. Mirzaei et al. [10] apply RL to
$\Delta$-tolling [11], which updates tolls in response to the difference
between the current travel time on each link and its free flow travel time.
Pandey et al. [5] apply DRL to congestion charging in a traffic model with a
mixture of managed and general lanes. Genser and Kouvelas [12] propose to use
multilayer neural networks to determine dynamic congestion tolls for inter-
area travel. Regardless, to our knowledge, DRL methods have not been utilized
to determine dynamic congestion tolls in general road networks with day-to-day
departure time and route choice problem.

##  III TRAFFIC MODEL

The definition of the traffic model considered in this study is as follows. We
consider a dynamic traffic model where departure time and route are determined
by day-to-day dynamics. The road network of the model has multiple Origin
Destination (OD) pairs connected by multiple routes. Each link may have a
bottleneck. Origins are residential areas, and destinations are central
business districts (CBD). The dynamic traffic model is an extension of the
bottleneck model [1].

Let $j$ be day, $t$ be time, $t^{\ast}$ be the travelers’ desired arrival time
to the destination. Let $\mu_{i}$ be the capacity of bottleneck $i$ per unit
time, $a_{j,i}(t)$ be the inflow rate at bottleneck $i$ on day $j$,
$N_{j,i}(t)$ be the number of vehicles in waiting queue on day $j$ and at
bottleneck $i$, $w_{j,i}(t)$ be waiting time on day $j$ and at bottleneck $i$,
and $\tau_{j,i}(t)$ be tolls on day $j$ and at bottleneck $i$. $t$ is time
when travelers leave bottleneck $i$ in them. In addition, let $B_{z}$ be a
given set of all bottlenecks on route $z$, and $Z_{\zeta}$ be a given set of
all routes connecting OD pair $\zeta$. Waiting queues are point queues with no
physical length.

The queue evolution at time $t$ is defined as Eq. (1).

$\frac{\mathrm{d}N_{j,i}(t)}{\mathrm{d}t}=\left\\{\begin{array}[]{ll}0\qquad(N_{j,i}(t)=0\
\mbox{and}\ a_{j,i}(t)<\mu_{i})\\\ a_{j,i}(t)-\mu_{i}\hskip
55.5pt(\mbox{otherwise})\end{array}\right.$ (1)

The waiting time in the queue at time $t$ is defined as Eq. (2).

$w_{j,i}(t)=\frac{N_{j,i}(t)}{\mu_{i}}$ (2)

The generalized cost of a traveler who arrives at the destination at time $t$
through route $z$ on day $j$ is defined as

$\begin{split}&c_{j,z}(t)=\left\\{\begin{array}[]{r}\sum\limits_{i\in
B_{z}}\tau_{j,i}(t_{\mathrm{l},z,i}(t))+\alpha\\{t-t_{\mathrm{d},z}(t)\\}+\beta\left(t^{\ast}-t\right)\\\
(t<t^{\ast})\\\ \\\ \sum\limits_{i\in
B_{z}}\tau_{j,i}(t_{\mathrm{l},z,i}(t))+\alpha\\{t-t_{\mathrm{d},z}(t)\\}+\gamma\left(t-t^{\ast}\right)\\\
(t\geq t^{\ast})\end{array}\right.\end{split}$ (3)

where $t_{\mathrm{d},z}(t)$ is a traveler’s departure time from an origin when
one arrives at a destination at time $t$ through route $z$,
$t_{\mathrm{l},z,i}(t)$ is a traveler’s departure time from bottleneck $i$
when one arrives at a destination at time $t$ through route $z$, $\alpha$ is
the value of time per unit time for a traveler, $\beta$ is the early arrival
penalty per unit time, and $\gamma$ is the late arrival penalty per unit time.

Multinominal logit model [13] is used for day-to-day dynamics. Travel cost is
defined with the weighted average learning operator [14] [15] as ”perceived
travel cost”. Note that its definition is extended in this study. The
probability of choosing departure time $t$ and route $z$ between OD pair
$\zeta$ is defined as

$P_{j,z}^{\zeta}(t)=\frac{\mathrm{exp}(-\vartheta\overline{c}_{j,z}(t))}{\sum\limits_{(z^{\prime},t^{\prime})}\mathrm{exp}(-\vartheta\overline{c}_{j,z^{\prime}}(t^{\prime}))}$
(4)

where $\vartheta$ is a parameter, and $\overline{c}_{j,z}(t)$ is the perceived
travel cost at day $j$, time $t$, and route $z$. Perceived travel cost is
defined as

$\begin{split}\overline{c}_{j,z}(t)=\frac{1}{\varsigma(\lambda)}\left(c_{j-1,z}(t)+\sum\limits_{i=2}^{T_{c}}\lambda^{i-1}c_{j-i,z}(t)\right)\\\
\forall z\in Z_{\zeta}\end{split}$ (5)
$\varsigma(\lambda)=\left\\{\begin{array}[]{ll}\sum_{i=1}^{T_{c}}\lambda^{i-1}&(0<\lambda<1)\\\
1&(\lambda=0)\end{array}\right.$ (6)

where $\lambda$ and $T_{\mathrm{c}}$ is a parameter, and
$1/\varsigma(\lambda)$ is a normalization factor.

Number of departing vehicles by departure time and route is defined as

$f_{j,z}(t)=M_{\zeta}\cdot P_{j,z}^{\zeta}(t)\qquad\forall z\in Z_{\zeta}$ (7)

where $f_{j,z}(t)$ is the number of departing vehicles by departure time and
route, and $M_{\zeta}$ is traffic demand between OD pair $\zeta$ per day.

Bounded rationality [13] is used with day-to-day dynamics. It means that the
departure time and route choice will not be changed on the next day when the
difference between the expected cost of the current choice and the minimum
expected cost of all other choices is less than or equal to a threshold
$\delta$. When a traveler who chooses route $z$ and departure time $t$ on day
$j-1$ does not change departure time and route choice on day $j$, the expected
cost $\widehat{c}_{j,z}(t)$ satisfies Eq. (8).

$\widehat{c}_{j,z}(t)-\widehat{c}_{j,z^{\prime}}(t^{\prime})\leq\delta\qquad\forall(z^{\prime},t^{\prime})\neq(z,t)$
(8)

##  IV CONGESTION PRICING METHOD

###  IV-A Overview

In this study, we propose ”Distributed Pricing Deep Deterministic Policy
Gradient (DP-DDPG)”, which is a DRL-based dynamic congestion pricing method.
Fig. 1 shows the outline of DRL.

Figure 1: The outline of DRL

In DP-DDPG, traffic data and toll update by bottleneck and time slot are
handled in DRL. In other words, a DRL agent performs ”distributed control”. It
significantly reduces the dimensionality of the states and actions.
Specifically, an DRL agent observes traffic data as a state at each pair of a
bottleneck and a time slot, and learns toll adjustment as an action at each
pair of a bottleneck and a time slot with reward.

In distributed control, cooperation among pairs of a bottleneck and a time
slot is required. In this study, we propose two methods for the cooperation:
spatially shared reward and temporally switching learning.

Deep Deterministic Policy Gradient (DDPG) [16] is used as the framework for
DRL. DDPG is one of DRL algorithms which deal with actions as continuous
values, uniquely decides an action for a state (deterministic policy), and is
based on actor-critic. Due to the space limitation, the detailed formulation
and algorithm of DDPG are omitted from this paper; interested readers may
refer the original paper [16].

The flow of DP-DDPG is shown as Fig. 2. Guo et al. [17] is used as its
reference.

Figure 2: The flow of DP-DDPG

###  IV-B Definition of state, action, and reward

A state is based from inflow rate to a bottleneck, waiting time and a toll.
The state at day $j$, time $t$, and bottleneck $i$ is defined as

$\begin{split}\bm{s}_{j,i}(t)=&\left(\frac{a_{j,i}(t)-\mu_{i}}{\mu_{i}},\frac{w_{j,i}(t)}{\sum\limits_{t\in
T}w_{0,i}(t)/\sum\limits_{t\in T}\varpi_{0,i,t}},\right.\\\
&\left.\quad\frac{\tau_{j,i}(t)-\bar{\tau}_{j,i}}{\sum\limits_{t\in
T}w_{0,i}(t)/\sum\limits_{t\in T}\varpi_{0,i,t}}\right)\end{split}$ (9)

$\varpi_{j,i,t}=\left\\{\begin{array}[]{ll}1&(w_{j,i}(t)>0)\\\
0&(w_{j,i}(t)=0)\end{array}\right.$ (10)

where $t$ is time when travelers leave bottleneck $i$, and $\bm{s}_{j,i}(t)$
is the state at day $j$, time $t$, and bottleneck $i$.

In order to input the state to neural networks, the scale of each element is
aligned by a capacity of a bottleneck and waiting time when no toll is charged
and day-to-day dynamics converge.

An action is a range of increase/decrease in a toll. The action at day $j$,
time $t$, and bottleneck $i$ is defined as

$b_{j,i}(t)=G\cdot\tanh y\quad(-G<b_{j,i}(t)<G)$ (11)

where $b_{j,i}(t)$ is the action at day $j$, time $t$, and bottleneck $i$, $y$
is an input to the final layer of an actor network, and $G$ is a parameter
which decides the range of an action.

The toll at day $j$, time $t$, and bottleneck $i$ is updated as Eq. (12).

$\tau_{j,i}(t)\leftarrow\tau_{j,i}(t)+b_{j,i}(t)$ (12)

A reward is constructed from the below two term.

  * •

a term based on waiting time by time slot and bottleneck

  * •

a term based on the average waiting time at all bottlenecks where tolls are
charged and all time slots (spatially shared reward)

The purpose of using spatially shared reward is cooperative decrease among all
bottlenecks where tolls are charged and all time slots. In a general road
network, fluctuation in waiting time at one bottleneck can affect that at
other bottlenecks. The term is expected to prevent from toll adjustment which
increases waiting time at other bottlenecks.

The reward at day $j$, time $t$, and bottleneck $i$ is defined as

$\begin{split}r_{j,i}(t)=&-\left\\{\frac{w_{j,i}(t)}{\sum\limits_{t\in
T}w_{0,i}(t)/\sum\limits_{t\in T}\varpi_{0,i,t}}\right.\\\
&\left.\qquad+\frac{1}{\\#I}\sum_{i\in
I}\frac{\bar{w}_{j,i}}{\sum\limits_{t\in T}w_{0,i}(t)/\sum\limits_{t\in
T}\varpi_{0,i,t}}\right\\}\end{split}$ (13)

where $r_{j,i}(t)$ is the reward at day $j$, time $t$, and bottleneck $i$, and
$I$ is a set of bottlenecks where tolls are charged.

###  IV-C Temporally switching learning

In this study, we propose ”temporally switching learning”. It means that an
agent does not perform DRL and fixes its action at zero in the pair of a time
slot $t^{\prime}$ and a bottleneck $i$ where a time slot moving average of
waiting time is less than a threshold $\delta$, and it is defined as

$\frac{1}{2n+1}\sum_{t^{\prime}=t-n}^{t+n}w_{j,i}(t^{\prime})<\mathrm{d}w$
(14)

where $n\in\mathbb{N}$ is a parameter, and $\mathrm{d}w$ is a constant.

The first purpose of temporally switching learning is improvement of the
efficiency of DRL. It narrows down the learning target of an agent to toll
adjustment in the pair of a time slot and a bottleneck where waiting time is
not nearly zero.

Its second purpose is prevention against excessive increase of tolls and
destabilization of learning in the pair of a time slot and a bottleneck where
DRL is turned off. In the pair of a time slot and a bottleneck where waiting
time keeps nearly zero, relatively high reward is given to an agent. This may
cause continuous increase of a toll and make learning unstable.

##  V NUMERICAL EXPERIMENTS

###  V-A Purpose and scenario

The performance of DP-DDPG is evaluated by various types of numerical
experiments. The summary of the experiment scenarios and their purposes are as
follows:

  * •

DP-DDPG is applied to Sioux Falls Network to evaluate the performance in the
condition where multiple OD pairs and multiple routes exist.

  * •

DP-DDPG is applied to a parallel bottleneck model, where route choice is
important, to compare it with other methods.

###  V-B Sioux Falls Network

We used Sioux Falls Network [18] as illustrated in Fig. 3.

Figure 3: Sioux Falls Network

Numbers of links with a bottleneck where a toll is charged and waiting time is
not zero are 29, 48, 53, and 58. Those with a bottleneck where a toll is not
charged and waiting time is not zero are 49, 52, and 61. $t^{\ast}$ is set to
be equal for all travelers.

DRL is performed in several patterns of values of learning rate of actor
networks and critic networks, and $G$. We decided that optimal parameters as
follows:

  * •

learning rate of actor networks: $10^{-3}$

  * •

learning rate of critic networks: $10^{-2}$

  * •

$G$: $1.5$

The parameters of traffic model are specified as follows: $\alpha=1.0$,
$\beta=0.45$, $\gamma=1.2$, $\vartheta=0.015$, $\lambda=0$, $\delta=1$,
$t\in\mathbb{N}$, $1\leq t\leq 250$, and $t^{\ast}=75$. During DRL, 40 days
are set as one cycle, and after 40 days have passed, DRL is returned to the
initial state. It is when ”$\tau_{j,i}(t)=0\hskip 10.0pt\forall(j,i,t)$ and
day-to-day dynamics converge”. Values of the parameters of actor networks and
critic networks are succeeded. One set of DRL consists of 15 cycles, and 10
sets are performed with the above values of learning rates and $G$.

Fig. 4 shows the transition of total travel time in Sioux Falls Network when
the trained agents in each set are applied to it under the same conditions.

Figure 4: Day-to-day dynamics of total travel time in Sioux Falls Network

On Fig. 4, initial total free travel time and initial total waiting time are
those when $\tau_{j,i}(t)=0\hskip 10.0pt\forall(j,i,t)$ and day-to-day
dynamics converge. According to Fig. 4, DP-DDPG stably reduces total travel
time in about 40 days.

Fig. 5 shows waiting time distribution and toll distribution at each
bottleneck.

Figure 5: Day-to-day dynamics of within-day waiting time (up) and toll (down)
distribution in each bottleneck

According to Fig. 5, tolls increase according to waiting time distributions.

###  V-C Parallel bottleneck model

We used a parallel bottleneck model with a single OD pair connected three path
and one bottleneck per path as illustrated in Fig. 6. It is the same model as
”multiple bottleneck model” [4]. $t^{\ast}$ is set to be equal for all
travelers.

Figure 6: Illustration of the parallel bottleneck model [4]

We compared the result of DP-DDPG with those of centralized DDPG, fully
distributed DDPG and Q-learning [4]. The outline of each method is as follows.
In centralized DDPG, a state is defined with traffic data at all bottlenecks
where tolls are charged, an action is defined with toll update there, and a
reward is based on waiting time there. Toll distributions are set as piecewise
linear functions. In other methods, a state is defined with traffic data at
each bottleneck where a toll is charged, and an action is defined with toll
update there. In DP-DDPG (our proposed method) and fully distributed DDPG, a
reward is based on waiting time at each bottleneck where a toll is charged. In
Q-learning, a reward is based on inflow rate at each bottleneck where a toll
is charged. In DP-DDPG and Q-learning, spatially shared reward is
incorporated.

Values of learning rate of actor networks and critic networks, and $G$ are as
follows:

  * •

learning rate of actor networks: $10^{-5}$

  * •

learning rate of critic networks: $10^{-4}$

  * •

$G$: $0.5$

The parameters of traffic model are specified as follows: $\alpha=1.0$,
$\beta=0.45$, $\gamma=1.2$, $\vartheta=0.05$, $\lambda=0$, $\delta=0.5$,
$t\in\mathbb{N}$, $1\leq t\leq 80$, and $t^{\ast}=30$.

During DRL by methods except Q-learning, 60 days are set as one cycle, and
after 60 days have passed, DRL is returned to the initial state. It is when
”$\tau_{j,i}(t)=0\hskip 10.0pt\forall(j,i,t)$ and day-to-day dynamics
converge”. Values of the parameters of actor networks and critic networks are
succeeded. One set of DRL consists of 20 cycles, and 20 sets are performed
with the above values of learning rates and $G$.

During RL by Q-learning, 1000 days are set as one cycle, and after 1000 days
have passed, RL is returned to the initial state. It is when
”$\tau_{j,i}(t)=0\hskip 10.0pt\forall(j,i,t)$ and day-to-day dynamics
converge”. The update result of parameter is succeeded. One set of RL consists
of 30 cycles, and three sets are performed. Note that replicator dynamics is
used as day-to-day dynamics in Q-learning.

Fig. 7 shows the transition of total waiting time in all methods.

Figure 7: Day-to-day dynamics of total waiting time in each method

According to Fig. 7, learned agents by centralized DDPG and those by fully
distributed DDPG could not decrease waiting time. In centralized DDPG, since
we used piecewise linear functions as toll distributions, they could not be
made close to the proper toll distributions, and day-to-day dynamics could not
become stable. In fully distributed DDPG, since cooperation among pairs of a
bottleneck and a time slot is not used, agents could not learn properly in the
condition where multiple bottlenecks exist. According to Fig. 7, learned
agents by Q-learning can decrease waiting time. It can be because spatially
shared reward was used. However, once waiting time became nearly zero, it
oscillated. It can be because actions are discrete values and temporally
switching learning is not used.

Note that waiting time reduction speed is not generally comparable because
replicator dynamics is used for the day-to-day dynamics and tolls are updated
every time day-to-day dynamics converge.

##  VI CONCLUSION

A dynamic congestion pricing method using deep reinforcement learning is
proposed. It is designed to eliminate traffic congestion based on observable
data in general large-scale road networks. One of the novel elements of the
proposed method is the distributed and cooperative learning scheme. It enables
a fast and computationally efficient learning in large-scale networks. The
numerical experiments using Sioux Falls Network showed that the method works
well thanks to the novel learning scheme.

## References

  * [1] W. Vickrey, “Congestion theory and transport investment,” _The American Economic Review_ , vol. 59, no. 2, pp. 251–260, 1969.
  * [2] H. Yang, Q. Meng, and D.-H. Lee, “Trial-and-error implementation of marginal-cost pricing on networks in the absence of demand functions,” _Transportation Research Part B: Methodological_ , vol. 38, pp. 477–493, 2004\.
  * [3] T. Seo and Y. Yin, “Optimal pricing for departure time choice problems with unknown preference and demand: Trial-and-error approach,” 2019.
  * [4] K. Sato, T. Seo, and T. Fuse, “A reinforcement learning-based dynamic congestion pricing method for the morning commute problems,” _Transportation Research Procedia_ , vol. 52, pp. 347–355, 2021, 23rd EURO Working Group on Transportation Meeting, EWGT 2020, 16-18 September 2020, Paphos, Cyprus.
  * [5] V. Pandey, E. Wang, and S. D. Boyles, “Deep reinforcement learning algorithm for dynamic pricing of express lanes with multiple access locations,” _Transportation Research Part C: Emerging Technologies_ , vol. 119, p. 102715, 2020.
  * [6] S. Wang, X. Xie, K. Huang, J. Zeng, and Z. Cai, “Deep reinforcement learning-based traffic signal control using high-resolution event-based data,” _Entropy_ , vol. 21, no. 8, 2019.
  * [7] Z. Li, H. Yu, G. Zhang, S. Dong, and C.-Z. Xu, “Network-wide traffic signal control optimization using a multi-agent deep reinforcement learning,” _Transportation Research Part C: Emerging Technologies_ , vol. 125, p. 103059, 2021.
  * [8] M. Li, Z. Li, C. Xu, and T. Liu, “Deep reinforcement learning-based vehicle driving strategy to reduce crash risks in traffic oscillations,” _Transportation Research Record_ , vol. 2674, no. 10, pp. 42–54, 2020. [Online]. Available: https://doi.org/10.1177/0361198120937976
  * [9] F. Zhu and S. V. Ukkusuri, “A reinforcement learning approach for distance-based dynamic tolling in the stochastic network environment,” _Journal of Advanced Transportation_ , vol. 49, pp. 247–266, 2015.
  * [10] H. Mirzaei, G. Sharon, S. Boyles, T. Givargis, and P. Stone, “Enhanced delta-tolling: Traffic optimization via policy gradient reinforcement learning,” in _2018 21st International Conference on Intelligent Transportation Systems (ITSC)_ , 2018, pp. 47–52.
  * [11] G. Sharon, M. W. Levin, J. P. Hanna, T. Rambha, S. D. Boyles, and P. Stone, “Network-wide adaptive tolling for connected and automated vehicles,” _Transportation Research Part C_ , vol. 84, pp. 142–157, 2017.
  * [12] A. Genser and A. Kouvelas, “Dynamic optimal congestion pricing in multi-region urban networks by application of a multi-layer-neural network,” _Transportation Research Part C: Emerging Technologies_ , vol. 134, p. 103485, 2022.
  * [13] Y. Yu, K. Han, and W. Ochieng, “Day-to-day dynamic traffic assignment with imperfect information, bounded rationality and information sharing,” _Transportation Research Part C: Emerging Technologies_ , vol. 114, pp. 59–83, 2020.
  * [14] E. Cascetta, “A stochastic process approach to the analysis of temporal dynamics in transportation networks,” _Transportation Research Part B: Methodological_ , vol. 23, no. 1, pp. 1–17, 1989.
  * [15] Y. Ouyang, “Pavement resurfacing planning for highway networks: Parametric policy iteration approach,” _Journal of Infrastructure Systems_ , vol. 13, no. 1, pp. 65–71, 2007.
  * [16] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning,” _arXiv preprint arXiv:1509.02971_ , 2015.
  * [17] S. Guo, X. Zhang, Y. Zheng, and Y. Du, “An autonomous path planning model for unmanned ships based on deep reinforcement learning,” _Sensors_ , vol. 20, no. 2, p. 426, 2020.
  * [18] E. Morlok, J. Schofer, W. Pierskalla, R. Marsten, S. Agarwal, J. Stoner, J. Edwards, L. LeBlanc, and D. Spacek, _Development and Application of a Highway Network Design Model, Volumes 1 and 2_ , ser. Final Report: FHWA Contract Number DOT-PH-11. Northwestern University, Evanston, 1973.


# 3D Neuron Morphology Analysis

Jiaxiang Jiang⋆ ⋆ Department of Electrical and Computer Engineering,
University of California, Santa Barbara
† Department of Molecular, Cell and Developmental Biology, University of
California, Santa Barbara  Michael Goebel⋆ ⋆ Department of Electrical and
Computer Engineering, University of California, Santa Barbara
† Department of Molecular, Cell and Developmental Biology, University of
California, Santa Barbara  Cezar Borba† ⋆ Department of Electrical and
Computer Engineering, University of California, Santa Barbara
† Department of Molecular, Cell and Developmental Biology, University of
California, Santa Barbara  William Smith† and B.S. Manjunath⋆ ⋆ Department of
Electrical and Computer Engineering, University of California, Santa Barbara
† Department of Molecular, Cell and Developmental Biology, University of
California, Santa Barbara

###### Abstract

We consider the problem of finding an accurate representation of neuron
shapes, extracting sub-cellular features, and classifying neurons based on
neuron shapes. In neuroscience research, the skeleton representation is often
used as a compact and abstract representation of neuron shapes. However,
existing methods are limited to getting and analyzing “curve”skeletons which
can only be applied for tubular shapes. This paper presents a 3D neuron
morphology analysis method for more general and complex neuron shapes. First,
we introduce the concept of skeleton mesh to represent general neuron shapes
and propose a novel method for computing mesh representations from 3D surface
point clouds. A skeleton graph is then obtained from skeleton mesh and is used
to extract sub-cellular features. Finally, an unsupervised learning method is
used to embed the skeleton graph for neuron classification. Extensive
experiment results are provided and demonstrate the robustness of our method
to analyze neuron morphology.

###### Index Terms:

3D Neuron Morphology, Skeleton Mesh, Graph, Sub-cellular features, Embedding,
Classification

##  I Introduction

The importance of neuronal morphology has been recognized from the early days
of neuroscience [1]. There are three obstacles in automatic neuron morphology
analysis. First, we need to have a good shape representation of each neuron.
Skeleton representations are widely used in neuroscience [2, 3, 4, 5] as they
provide a compact and abstract shape representation. Mathematically,
skeletonization or medial axis transform (MAT) has a rigorous definition for
arbitrary shapes. The skeleton of a shape is defined as a collection of
interior points that have at least two closest points on the surface of the
shape. We refer to those interior points as skeleton points and each skeleton
point is associated with a radius. Figure 1 shows an example of MAT. However,
in reality, it is not an easy task to get skeleton representation directly
from images. Most automatic or manual segmentation methods output a cloud of
surface points. Thus, we need to compute the 3D neuron skeleton from 3D
surface point clouds. The skeleton representation further enables computing
sub-cellular features such as length and number of branches of neurons as well
as classification of neurons.

Figure 1: Illustration of MAT by using a 2D shape example.

The main contribution of this paper is a robust and efficient method for
computing a skeleton representation from a set of 3D surface points. This 3D
skeleton representation can be used for a quantitative analysis of neuronal
cell structures, including sub-cellular feature calculations and for neuron
type classification based on 3D shapes.

There is an extensive literature on neuron skeleton extraction [3, 6, 7]. In
[3], the skeleton representation is computed from a 3D mesh by using a
traditional morphological thinning algorithm [8]. This method has two main
drawbacks. First, the thinning algorithm is sensitive to noise of 3D mesh.
Second, in reality, we usually get discrete 3D surface points of neurons from
the segmentation step, and constructing 3D mesh from those discrete 3D surface
points will introduce additional noise. To make the skeleton extraction model
more robust, [7] proposes to use deep learning network to learn skeleton
representation. The main idea of the paper is to use the deep learning network
to predict skeletons from features of multiple spatial scale layers. This
model still takes a continuous surface as input, as opposed to discrete
surface points. Further, this is a supervised method and it requires training
samples. In [6], they propose extracting skeleton representations directly
from discrete surface points by using a 3D discrete distance transform.
However, this only works well for curve skeletons and only tubular structures
have curve skeletons. General 3D shapes will result in surface skeletons as
shown in Fig.2.

Figure 2: Visualization of a 3D ellipsoid shape and its surface skeleton from
two points of view. Yellow triangle mesh represents object surface. Black
contour represents the outline of the skeleton surface. Magenta and Cyan line
segments represent two closest surface points from the skeleton point. Two
colors are used to differentiate different directions.

In practice, skeleton mesh is used to represent the surface skeleton. There is
no existing method to extract skeleton mesh from surface point clouds for
neuron morphology analysis.

There are also methods to analyze neuron shapes using the skeleton
representation. For skeleton classification task, [3] proposes to compare
similarity of skeletons by using local skeleton features. It breaks a neuron
skeleton into short segments and characterize segments by location and
direction of segments. However, this method only works well for skeleton
curves but not surface skeletons. In [9], they only consider features from
skeleton points for the classification and it does not fully utilize the
skeleton mesh information.

To analyze general neuron shapes, this paper presents a robust 3D neuron
morphology analysis framework based on the surface skeleton representation of
neurons. In [10], the authors propose an unsupervised deep learning skeleton
mesh extraction method. However, this method does not work well when neurons
have concave shapes. Our skeleton mesh extraction method is built upon [10],
and by using estimated surface norm of point clouds as part of the
optimization function, we address this drawback. Next the skeleton mesh is
converted to an undirected graph called skeleton graph. Inspired by [11], we
embed the skeleton graph by maximizing mutual information, and then classify
neurons based on the embedding of each skeleton. To compute cellular/sub-
cellular features of neurons from the skeleton representation, we also utilize
the skeleton graph. A simple but effective recursive algorithm is proposed to
get number of branches and length of neurons.

We apply our neuron morphology analysis method to classify Ciona neurons. The
Ciona sea squirt is one of the widely studied tunicates in neuroscience [12].
Its brain is closely related to vertebrates with a much simpler neuronal
structure. In a single Ciona larva, it has only about 187 neurons with about
6600 synapses [12]. Studying the Ciona brain in depth can reveal the general
principles behind the mechanism of how vertebrate brains work [13]. We also
present our results on the NeuroMorpho[14] dataset. In summary, the main
contribution of our paper include

  * •

A robust and efficient skeleton mesh extraction method with novel cost
function by using properties of MAT. To the best of our knowledge, this is the
first one to use skeleton mesh instead of the curve skeleton to analyze
neuronal shapes

  * •

A 3D Ciona neuron dataset that can be used for neuron morphology analysis.

##  II Method

Figure 3: Overview of our proposed neuron morphology analysis pipeline. Given
a surface point cloud as input, we extract the skeleton mesh. The skeleton
mesh includes skeleton points with their radii as well as the connection of
skeleton points. Then we construct the skeleton graph. Each node in a skeleton
graph represents a skeleton point, and edge in the graph represents the
connection between skeleton points. Next, we propose a graph analysis method
to get length and number of branches of neurons based on the skeleton graph.
We also use the skeleton graph for classification task by embedding it into a
fixed length vector.

Fig.3 illustrates our overall neuron morphology analysis method. Given a set
of surface point clouds as input, we introduce an unsupervised deep learning
method to get the skeleton mesh representation of each neuron. This is
achieved by using the properties of the traditional medial axis transform
(MAT). The skeleton mesh representation includes skeleton points with radii as
well as the connection of those skeleton points as shown in Fig.3. Second, the
skeleton representation of each neuron is transformed into a skeleton graph.
Each node in the skeleton graph represents a skeleton point. If there is an
edge between two nodes, it means those two skeleton points are connected. The
weight of the edge represents distance between the two skeleton points. Radii
as well as the location of each skeleton point are attributes of each node.
Next, length and number of branches of neurons are computed from the skeleton
graph. To compare different shapes of neurons, a graph level representation
learning method is used to embed the skeleton graph. The representation
learning method is an unsupervised method that maximizes mutual information of
the skeleton graph.

###  II-A Skeleton Mesh from 3D Surface Point Cloud

Figure 4: Overview of neuron skeleton representation method. Given 3D surface
point cloud as input, PointNet++[15] is used to extract features of the input
point cloud. Then a geometric transformation learned by MLP will predict the
skeleton points location with their radii. After skeleton points prediction,
two simple priors are used to initially connect some skeleton points, and a
graph auto-encoder is used to predict all links that connect skeleton points.

Our unsupervised 3D neuron skeleton extraction method is built upon the method
in [10] as illustrated in Fig.4. Given a 3D surface point cloud as input,
PointNet++ [15] is used as the encoder to obtain the sampled surface points
with features. Next, a multi-layer-perceptron (MLP) is used to learn the
geometric transformation to predict the skeleton points with their radii using
linear combination of the MLP input points. Compared to [10], we propose to
use properties of skeleton points as the prior knowledge to make the geometric
transformation learning process more robust to general shapes. After getting
skeleton points with radii, a graph auto-encoder is used to predict links
between skeleton points.

####  II-A1 skeleton points prediction

Mathematically, given a set of 3D surface points
$\textbf{P}\in\mathbb{R}^{M\times 3}$ where $M$ is a number of points, we want
to predict $N$ skeleton spheres $\textbf{s}_{i}=<\textbf{c}_{i},r_{i}>$ where
$\textbf{c}_{i}$ is the center of each sphere and $r_{i}$ is the radius of
each sphere.

As illustrated in Fig.4, we first use PointNet++[15] as the encoder to obtain
a set of sampled input points
$\textbf{P}^{\prime}\in\mathbb{R}^{M^{\prime}\times 3}$ and their associated
contextual features $\textbf{F}\in\mathbb{R}^{M^{\prime}\times D}$.
$M^{\prime}(M^{\prime}<M)$ is the number of feature points after PointNet++
and $D$ is the feature dimension of each feature point. PointNet++ groups
points and extract point features hierarchically. It contains a number of set
abstraction levels. For each set abstraction level, there are three layers:
Sampling layer, Grouping layer, and PointNet layer. For the first set
abstraction level, the input is P, a set of $M$ number of 3D surface points.
Next, Sampling layer is applied. The iterative farthest point sampling (FPS)
[15] is used to get $M^{\prime}$ sampled points. In the grouping layer, those
$M^{\prime}$ sampled points are used as the centroid points. Then for each
centroid, all $M$ points within a radius are viewed as neighbor points and are
grouped into that local region. Therefore, each centroid has $K$ neighbors and
$K$ can vary for different groups. After Sampling and Grouping layers,
PointNet[16] layer is used to extract features for each local region. The
sampling layer, grouping layer, and PointNet layer consists one set
abstraction level, and we stack such abstraction levels to form a hierarchical
architecture to get features at different spatial scales. Next, multi-scale
grouping is applied to concatenate the features from different spatial scales.

A Multi-Layer Perceptron (MLP) is used to estimate the geometrical
transformation to get a set of skeleton spheres’ center points C,
$\\{\textbf{c}_{1},\textbf{c}_{2},...,\textbf{c}_{N}\\}$. The geometric
transformation we use is a convex combination of input points
$\textbf{P}^{\prime}$. MLP with softmax function is used to estimate the
weight $\textbf{W}\in\mathbb{R}^{M^{\prime}\times N}$ of the convex
combination in eq. 1.

$\textbf{C}=\textbf{W}^{T}\textbf{P}^{\prime}\qquad\text{subject
to}\qquad\forall j\in\\{1,...,N\\}\quad\sum_{i=1}^{M^{\prime}}W_{i,j}=1$ (1)

As shown in [10], the same weight matrix W can be used to compute
$r(\textbf{c})\in\textbf{R}$ using eq.2

$\textbf{R}=\textbf{W}^{T}\textbf{D}$ (2)

where $D\in\mathbb{R}^{M^{\prime}\times 1}$ is a vector of
$d(\textbf{p}^{\prime},\textbf{C})$. $d(\textbf{p}^{\prime},\textbf{C})$ is
the closest distance of one surface point $\textbf{p}^{\prime}$ to all
skeleton points C and is defined as
$d(\textbf{p}^{\prime},\textbf{C})=\min_{\textbf{c}\in\textbf{C}}||\textbf{p}^{\prime}-\textbf{c}||_{2}$

A set of loss functions are defined in [10] to train the MLP. The loss
function includes sampling loss $L_{s}$, point-to-sphere loss $L_{p}$, and
radius regularizer loss $L_{r}$. The first two losses are based on the
recoverability of skeleton representation. The last loss term is to encourage
larger radii to avoid instability induced by surface noise.

For the sampling loss $L_{s}$, we sample points on the surface of each
skeleton sphere and measure the Chamfer Distance (CD) between the set of
sampled sphere points T and the set of sampled surface points from PointNet++
$\textbf{P}^{\prime}$ as in eq. 3:

$L_{s}=\sum_{\textbf{p}^{\prime}\in\textbf{P}^{\prime}}\min_{\textbf{t}\in\textbf{T}}||\textbf{p}^{\prime}-\textbf{t}||_{2}+\sum_{\textbf{t}\in\textbf{T}}\min_{\textbf{p}^{\prime}\in\textbf{P}^{\prime}}||\textbf{t}-\textbf{p}^{\prime}||_{2}$
(3)

Point-to-sphere loss $L_{p}$ measures the reconstruction error by explicitly
optimizing the coordinate of skeleton points and their radii:

$L_{p}=\sum_{\textbf{p}^{\prime}\in\textbf{P}^{\prime}}(\min_{\textbf{c}\in\textbf{C}}||\textbf{p}^{\prime}-\textbf{c}||_{2}-r(\textbf{c}_{p^{\prime}}^{min}))+\sum_{\textbf{c}\in\textbf{C}}(\min_{\textbf{p}^{\prime}\in\textbf{P}^{\prime}}||\textbf{c}-\textbf{p}^{\prime}||_{2}-r(\textbf{c}))$
(4)

where C is a set of predicted skeleton points, $r(c)$ is a radius of skeleton
point c, and $\textbf{c}_{p^{\prime}}^{min}$ is the cloest skeleton points to
a point $\textbf{p}^{\prime}$.

Radius regularizer loss $L_{r}$ is defined in eq.5 where $r(\textbf{c})$ is a
radius of the skeleton point c. By minimizing this loss, it encourages large
radii of skeleton points to make the skeleton points prediction more stable.

$L_{r}=-\sum_{\textbf{c}\in\textbf{C}}r(\textbf{c})$ (5)

However, based on above three losses, predicted skeleton points can be outside
of a shape if the shape is concave. Therefore, we introduce the skeleton-to-
surface norm loss $L_{n}$ to encourage the skeleton points to be inside the
shape. $L_{n}$ is a term to measure the reconstruction error by utilizing the
property of spoke direction in MAT. Fig.5 illustrates a spoke of a skeleton
point in MAT. The length of a spoke of the skeleton point c is
$r(\textbf{c})$. This is also one of our main contribution compared to [10]
for skeleton points prediction. In theory, the direction of the spoke should
be perpendicular to the object surface at the surface point [9]. Also the
spoke direction should be pointing outside of a shape. To capture this
property, we define $L_{n}$:

$L_{n}=\sum_{\textbf{c}\in\textbf{C}}(1-\textbf{n}_{\textbf{p}_{c}^{{}^{\prime}min}}\cdot\frac{\textbf{p}_{c}^{{}^{\prime}min}-\textbf{c}}{||\textbf{p}_{c}^{{}^{\prime}min}-\textbf{c}||_{2}})+\sum_{\textbf{p}^{\prime}\in\textbf{P}^{\prime}}(1-\textbf{n}_{p^{\prime}}\cdot\frac{\textbf{p}^{\prime}-\textbf{c}_{p^{\prime}}^{min}}{||\textbf{p}^{\prime}-\textbf{c}_{p^{\prime}}^{min}||_{2}})$
(6)

$\textbf{p}_{c}^{{}^{\prime}min}$ is the closest surface point to the skeleton
point $c$ and $\textbf{n}_{\textbf{p}_{c}^{{}^{\prime}min}}$ is the estimated
surface norm of the surface point. The “$\cdot$” denotes the dot product
between two vectors. To estimate the norm of each point in the 3D surface
points, the adjacent points are found first and then principal axis of the
adjacent points using covariance analysis are calculated. More details of the
norm estimation of each surface point are described in [17]. $L_{n}$
encourages the skeleton points to be located within a shape even the shape is
concave.

Figure 5: Spoke is a vector connecting a skeleton point and that skeleton
point’s one of two closest surface points. The vector points from the skeleton
point to the surface point. Green dashed lines represent the implicit surface
of an object, blue dot is one skeleton point, orange dots represent surface
points, and the arrow represents a spoke. Spoke is perpendicular to the object
surface at the surface point.

####  II-A2 Links prediction

After predicting skeleton points, our target is to predict links to connect
skeleton points to form the skeleton mesh. In theory, for any pair of skeleton
points, if all points that are on the line connecting the two skeleton points
are also on the skeleton surface, there should be links between those two
points. We adapt the graph auto encoder (GAE) as used in [10] to predict links
between skeleton points. GAE takes input the initialized adjacency matrix
$A_{ini}$ of the skeleton mesh graph $G_{mesh}$ and the skeleton points’
features. The skeleton points’ features is concatenation of C, R, and
$\textbf{W}^{T}\textbf{F}$. C are coordinates of skeleton points, R are radii
of skeleton points, W is the learned weights from the MLP, and F is the
contextual features of the surface points from PointNet++. GAE outputs the
estimated adjacency matrix $\hat{A}$ of $G_{mesh}$. The loss function is a
Masked Balanced Cross-Entropy (MBCE) loss as proposed in [18].

###  II-B Sub-cellular feature extraction from the skeleton graph

Skeleton model can be represented as the skeleton graph $G(V,E)$ where $V$
represents all skeleton points and $E$ represents connection between skeleton
points. The weight of the edge represents distance between skeleton points.

#### Neuron length computation from skeleton graph

We formulate the neuron length computation problem as finding the longest
simple path in the skeleton graph $G$. A simple path in the graph is a path
that does not have repeat nodes. In $G$, the length of the path is the sum of
all edges’ weights along the path. Note that our skeleton graph can have
loops. To solve this problem, for each node, we find the longest simple path
from that node and denote as $path(i)$ for node $i$. To avoid getting stuck
during the loop, we mark any node when we visits as shown in the algorithm
below. Next, we find $i^{*}$ that maximizes $path$. We use the following
recurrent algorithm to find the longest simple path from one node in the
skeleton graph.

Algorithm 1 Find the longest simple path from node i

D[j] represents the longest path from j to i. Initially, D[j]=0 for all j

function LongestPath(i, currLength)

if Node i is visited then

return

end if

change i status to visited

if D[i] <currLength then

D[i]=currLength

end if

for all nodes j that is connected to i do

LongestPath(j, currLength+edgeweight[i,j])

end for

change i status to not visited

end function

#### Neuron branch calculation

After finding the longest simple path, we are able to identify a set of nodes
on that path. Those nodes are possible branching nodes. We name a set
containing all possible branching nodes as B For each node $i\in B$, we find
the longest simple path from that node $i$ which does not contain any other
nodes in $B$. Therefore, the branch is identified as the longest simple path.

###  II-C Skeleton Model Comparison

Figure 6: We use two example skeleton graphs (blue and orange) to demonstrate
how we embed the skeleton graph. Each node of a skeleton graph is encoded into
a feature vector by using graph convolution layers. A fixed length graph level
feature vector (global representation) is obtained by graph-level pooling
operation of each node feature vector. The discriminator takes inputs both
global representation and patch representation to decide whether they are from
the same skeleton graph. In this toy example, there will be 14 global-patch
pairs.

We cluster neuron morphology by comparing different skeleton graphs.
Specifically, we embed the skeleton graphs and then cluster neurons based on
their embeddings. We embed the skeleton graph based on InfoGraph [11] as
illustrated in Fig. 6. The embedding process is in an unsupervised manner.

First, graph convolutional layers are used to generate node features which we
refer to as patch representation $\textbf{h}_{i}^{j}$ ($i$ is the skeleton
graph index and $j$ is the node index of the skeleton graph $i$). Then graph-
level pooling is used on all patch representations to get the graph level
representation (global representation) $\textbf{H}_{i}$. The mutual
information (MI) estimator on global-patch pairs over the given graph dataset
G is defined as:

$MI=\sum_{i\in K}\frac{1}{K}\sum_{j\in
G_{i}}I(\textbf{h}_{i}^{j},\textbf{H}_{i})$ (7)

where $K$ is the total number of graphs in the dataset and
$G_{i}\in\textbf{G}$.

MI is the mutual information estimator modeled by the discriminator $T$. The
Jensen-Shannon MI estimator proposed in [11] is

$I(\textbf{h}_{i}^{j},\textbf{H}_{i})=\mathbb{E}[-sp(-T(\textbf{h}_{i}^{j},H_{i})]-\mathbb{E}[-sp(-T(\textbf{h}_{i^{\prime}}^{j},H_{i})]$
(8)

where $\mathbb{E}$ is the expectation (here it is just average operation) and
$sp(z)=log(1+e^{z})$. $i$ and $i^{\prime}$ denote two graph samples from the
dataset G. The discriminator $T$ estimates global-patch representation pairs
by passing two representations to different non-linear transformations and
then takes the dot product of the two transformed representations. Both non-
linear transformations consist 3 linear layers with ReLU activation functions.
Therefore, the discriminator will output a score between $[0,\infty)$ to
represent whether the input patch (node) is from the input graph. If the input
global/patch pairs are from the same graph, we refer to them as positive
samples, otherwise negative samples. We randomly sample pairs as input to the
discriminator.

##  III Dataset

###  III-A Ciona Neuron EM Dataset

The first dataset (Dataset 1) contains two Ciona larva 3D TEM images [12]. The
section thickness for TEM images varies between 35 nm and 100 nm. For each
section, xy resolution is 3.85$\times$3.85 nm. Animal 1 contains 7671 sections
and animal 2 contains about 8000 sections. In each Ciona larva, 187 neurons
are annotated. Those 187 neurons can be grouped into 31 types. For animal 2,
Ciona neuron skeletons are traced using TrackEM2 [19], an ImageJ [20] plugin.
This dataset is summarized in Table I, and we refer to [12] for more details.

TABLE I: Details of Ciona Dataset. It contains two Ciona animals, one with surface point cloud annotated and one with skeleton annotated.  Animal | xy resolution (nm) | section thickness (nm) | number of sections | annotations
---|---|---|---|---
Animal 1 |  3.85$\times$3.85 | 35-100 | 7671 | 3D surface point cloud of neurons are provided
Animal 2 |  3.85$\times$3.85 | 35-100 | 6928 | 3D neuron skeletons without skeleton points’ radii

###  III-B C.elegans Neuron Dataset from NeuroMorpho

NeuronMorpho [14] is a publicly available dataset that is used for neuron
morphology research. It has dozens of different animals’ neurons. So far, it
is the largest neuron skeletons dataset with associated metadata. In this
paper, we take a subset of C.elegans dataset (Dataset 2) from the whole
NeuroMorpho dataset to verify our sub-cellular feature extraction method.
Dataset 2 consists of 299 neuron skeletons (with radii) and it is classified
into 10 different types. Each neuron with detailed metadata information such
as number of branches and length of neuron.

##  IV Results

###  IV-A Skeleton Model from 3D Surface Point Cloud

We apply our method on animal 1 neurons from Dataset 1 for the purpose of
building a shape model to analyze neuron morphology. To get the fixed number
of 3D surface input points, we use the sampling strategy described in [21].
The main idea of the sampling strategy is to give each point a weight based on
its distance to neighbor points. Then we sample points based on the weights
until we reach the number of desired points. Details of defining the neighbor
points and computing the weights are described in [21].

\begin{overpic}[width=86.72267pt]{images/skeleton_mesh/2_raw.png}\put(0.0,80.0){A}\end{overpic}

\begin{overpic}[width=86.72267pt]{images/skeleton_mesh/2_DPC.png}\put(0.0,80.0){B}\end{overpic}

\begin{overpic}[width=86.72267pt]{images/skeleton_mesh/2_point2skeleton.png}\put(0.0,80.0){C}\end{overpic}

\begin{overpic}[width=86.72267pt]{images/skeleton_mesh/2_our_method.png}\put(0.0,80.0){D}\end{overpic}

\begin{overpic}[width=86.72267pt]{images/skeleton_mesh/pr13_raw.png}\put(0.0,30.0){A}\end{overpic}

\begin{overpic}[width=86.72267pt]{images/skeleton_mesh/pr13_DPC.png}\put(0.0,30.0){B}\end{overpic}

\begin{overpic}[width=86.72267pt]{images/skeleton_mesh/pr13_point2skeleton.png}\put(0.0,30.0){C}\end{overpic}

\begin{overpic}[width=86.72267pt]{images/skeleton_mesh/pr13_our_method.png}\put(0.0,30.0){D}\end{overpic}

\begin{overpic}[width=86.72267pt]{images/skeleton_mesh/21_raw.png}\put(0.0,80.0){A}\end{overpic}

\begin{overpic}[width=86.72267pt]{images/skeleton_mesh/21_DPC.png}\put(0.0,80.0){B}\end{overpic}

\begin{overpic}[width=86.72267pt]{images/skeleton_mesh/21_point2skeleton.png}\put(0.0,80.0){C}\end{overpic}

\begin{overpic}[width=86.72267pt]{images/skeleton_mesh/21_our_method.png}\put(0.0,80.0){D}\end{overpic}

\begin{overpic}[width=86.72267pt]{images/skeleton_mesh/AMG3_raw.png}\put(0.0,40.0){A}\end{overpic}

\begin{overpic}[width=86.72267pt]{images/skeleton_mesh/AMG3_DPC.png}\put(0.0,40.0){B}\end{overpic}

\begin{overpic}[width=86.72267pt]{images/skeleton_mesh/AMG3_point2skeleton.png}\put(0.0,40.0){C}\end{overpic}

\begin{overpic}[width=86.72267pt]{images/skeleton_mesh/AMG3_our_method.png}\put(0.0,40.0){D}\end{overpic}

\begin{overpic}[width=86.72267pt]{images/skeleton_mesh/ant1_raw.png}\put(0.0,30.0){A}\end{overpic}

\begin{overpic}[width=86.72267pt]{images/skeleton_mesh/ant1_DPC.png}\put(0.0,30.0){B}\end{overpic}

\begin{overpic}[width=86.72267pt]{images/skeleton_mesh/ant1_point2skeleton.png}\put(0.0,30.0){C}\end{overpic}

\begin{overpic}[width=86.72267pt]{images/skeleton_mesh/ant1_our_method.png}\put(0.0,30.0){D}\end{overpic}

Figure 7: The figure shows skeleton extraction results from different methods.
A. Input 3D surface points; B. sampled skeleton points from surface points
using DPC [22]; C. skeleton mesh from surface points using Point2Skeleton
[10]; D. skeleton mesh from our method with surface norm cost function.

To evaluate our skeleton extraction method on Dataset 1, we carefully repair
surface mesh using screened poisson surface reconstruction method [23] with
spherical harmonics to smooth the surface. Fig.7 shows the qualitative
comparison of our methods and other state-of-the-art methods [22, 10]. Our
method has better visual results. [22] can generate unstructured skeleton
points but it lacks topological constraint. We sample the number of points
using [21] to be the same with the other methods for a fair comparison. It
performs well when neuron has tube like structure but it is not good when
neuron has a more circular shape. Compared with [10], our method can capture
more detailed structures which are important for sub-cellular feature
extraction, such as branches. For a quantitative evaluation of our method on
Dataset 1, we compute the strictly defined MAT and use the handcrafted method
in [10] to remove insignificant spikes to get the simplified MAT. We sample
points on the simplified MAT as the ground truth skeleton points. Then we
compute Chamfer Distance (CD) and Hausdorff distance (HD) between computed
skeleton points and ground truth skeleton points. We refer to them as CD-skel
and HD-skel, respectively. To compute CD-skel and HD-skel, we shift and
rescale each skeleton so that the skeleton center is located at (0,0,0) and
their x,y,z coordinates are all between -1 and 1 for every skeleton. We also
measure the difference between the shapes reconstructed from the skeletons and
ground truth surface points using CD and HD. We refer to them as CD-recon and
HD-recon, respectively. Similarly, we also shift and rescale the ground truth
points so that each neuron is centered at (0,0,0) and each neuron’s surface
coordinates are between -1 and 1. Other than those four aforementioned
evaluation metrics, we also use the reconstructed neuron volume difference as
the evaluation metric, considering neuron volume as one of the important
property of a neuron. We denote it as $vol-pct$. Mathematically, it is defined
as $vol-pct=\frac{|vrecon-vg|}{vg}$ where $vrecon$ is the volume of the
reconstructed neurons from the skeleton model, and $vg$ is the ground truth
volume. Table II gives the detailed results of different methods. It shows our
method has the best performance compared to other methods on Dataset 1 in
terms of all of the above evaluation metrics.

TABLE II: Quantitative Comparison with state-of-the-art skeleton model extraction method on Dataset 1 str | CD-recon | HD-recon | CD-skel | HD-skel |  vol-pct ($\%$)
---|---|---|---|---|---
DPC [22] | 0.102 | 0.298 | 0.303 | 0.311 | 10.1
Point2Skeleton [10] | 0.081 | 0.207 | 0.155 | 0.191 | 8.2
Our method | 0.067 | 0.183 | 0.090 | 0.185 | 5.6

###  IV-B Sub-cellular feature extraction from skeleton model

We apply our sub-cellular feature extraction method on Dataset 1 and Dataset
2. We define the length difference percentage (len-pct) and number of branches
difference percentage (num-pct) to measure the neuron length error and number
of branches error of the computation methods. We compare our sub-cellular
feature extraction method with the state-of-the-art sub-cellular feature
extraction method proposed in [3]. Table III gives details of sub-cellular
feature computation results. len-pct and num-pct for Dataset 1 is only for
both animals. Our method provides the better sub-cellular feature extraction
results in most cases and the percentage error is no more than 8 percent.
Also, as we see, our method is comparable to [3] on ”curve” skeletons but has
better performance on the skeleton meshes (”curve” skeleton is a special case
of the skeleton mesh).

TABLE III: Sub-cellular feature evaluation results |  len-pct on Dataset 1 ($\%$) |  num-pct on Dataset 1 ($\%$) |  len-pct on Dataset 2 ($\%$) |  num-pct on Dataset 2 ($\%$)
---|---|---|---|---
NAVIS[3] | 7.1 | 10.5 | 3.2 | 6.3
Our method | 5.5 | 7.6 | 3.8 | 5.8

For Dataset 1, we also analyze the relationships between length and branches
of neurons computed from our method. For animal 1 in Dataset 1, we first use
our skeleton extraction method to get the skeleton representation of each
neuron. Next, we use our sub-cellular feature extraction method to get length
of the neuron and number of branches of the neuron. Fig. 8 shows the
relationships between length and branches of neurons. Blue dots represent
neurons from animal 1 and red dots represent neurons from animal 2. Overall,
longer neurons tend to have more branches. This relationship between neuron
length and number of branches is consistent between animals 1 and 2.

Figure 8: Relationships between length and number of branches of neurons using
two animals of Dataset 1. Blue dots represent neurons from animal 1 and red
dots represent neurons from animal 2.

###  IV-C Skeleton Model Comparison

We apply our skeleton model comparison method on Dataset 1 for the purpose of
analyzing Ciona neuron morphology of different neuron types. Given the
skeleton graph, we embed it into a 100-dimension vector. Then we use K-means++
to cluster vector representations of skeleton graphs. For K-means++, the
number of clusters is set to be the same as number of neuron classes. After
K-means++, the cluster label is given by using the majority vote of neuron
types within the cluster. We use animal 1 neurons to train the K-means++ model
and get the cluster (neuron type) centers. Next, we use animal 2 as the test
set. For each neuron in the test set, we assign the label based on its closest
cluster center. The distance metric we use is the euclidean distance in the
vector embedding space. Table IV shows the comparison of clustering
(classification) accuracy on both training and test sets using different
neuron classification method. The neuron classification methods include graph
spectrum method, graph2vec [24], s-rep [9] and our graph level representation
method. The graph spectrum method uses the eigen values of the graph’s
adjacency matrix to form the vector representation. Similar to our method,
graph2vec method is another way to convert the skeleton graph to the graph
level vector representation. For the s-rep method, it uses the skeleton
points’ features such as spoke direction, spoke length, and skeleton points’
locations to classify neurons. From Table IV, we observe that grouping neurons
by our graph embedding provide the best classification results on both train
and test sets. It shows that neuron types are closely related to its
morphology. Also, our method is a better way to represent skeleton graphs in
terms of clustering accuracy.

TABLE IV: Neuron Classification Results | Graph Spectrum |  Graph2vec[24] |  S-rep[9] | Our method
---|---|---|---|---
Train | 0.691 | 0.767 | 0.791 | 0.893
Test | 0.632 | 0.718 | 0.773 | 0.871

Based on previous observations, we do further morphology analysis based on our
graph level representation results. After we get the vector representation of
each graph, we compute euclidean distance between each pair of vectors. Then
we compute the inter class and intra class distance based on pairwise neuron
distance as Fig 9 shows. Diagonal entries tend to be smaller than other
values, confirming a strong correlation between structure and function. More
specifically, neurons within a neuron type tend to have a smaller
morphological distance than neurons between different groups. Also, two
animals inter and intra distance look very similar.

Based on this inter and intra class distance, we do hierarchical clustering as
shown in Fig 10. The hierarchical clustering results show that BVIN and pr-BTN
RN have larger morphology distances from other neuron types. The BVIN neurons
are a broad group of intrinsic interneurons located in the brain vesicle of
Ciona. The main role of this group is to connect the sensory neurons to other
groups within the brain vesicle. The BVIN neurons have partial
subclassification based on sensory input [12]. Receiving specific sensory
information is an indication of functional role, therefore, the BVIN can be
further subdivided into different groups based on the sensory group(s) from
which they receive input. Using the sensory input as criteria, the entire
group was split up into four groups: prIN if receiving photoreceptor input,
antIN if receiving antenna cell input, pr-ant IN if receiving from both, and
BVIN if not receiving from either. The pr-BTN RN only have two neurons and
their functions are mostly unknown. According to the connectome [12], they
receive input from both the photoreceptors and the BTN neurons (neurons
involved in processing touch stimuli in the tail), so it’s possible they play
a role in integrating the two inputs. Any functional differences that may
exist between the two are currently unknown, however, the hierarchical
clustering suggests that this may be the case.

\begin{overpic}[width=260.17464pt]{images/pair_wise_log10_animal1.png}\put(5.0,4.0){A}\end{overpic}

\begin{overpic}[width=260.17464pt]{images/pair_wise_log10_animal2.png}\put(5.0,4.0){B}\end{overpic}

Figure 9: Inter and intra class neuron morphology distance on animal 1 (A) and
animal 2 (B) . Neuron morphology distance is computed by using euclidean
distance between our graph level representation of the skeleton graph.

\begin{overpic}[width=260.17464pt]{images/distance_cluster_animal1.png}\put(5.0,4.0){A}\end{overpic}

\begin{overpic}[width=260.17464pt]{images/distance_cluster_animal2.png}\put(5.0,4.0){B}\end{overpic}

Figure 10: Hierarchical clustering of neurons of animal 1 (A) and animal 2
(B).

##  V Conclusion

In this paper, we propose a novel neuron morphology analysis pipeline. It
mainly includes three parts. First, we propose a robust shapre representation
using skeleton mesh. Next, we compute sub-cellular features from the skeleton
mesh. Finally, we compare different neuron shapes using skeleton mesh. To the
best of knowledge, this is the first time that such an approach is used to
represent and classify neuronal shapes. The introduction of the estimated
surface norm penalty results in a robust mesh representation that achieves the
state-of-the-art performance using well defined metrics. Based on skeleton
graph, we formulate sub-cellular feature computation problem as a longest
simple path problem that can be easily computed. To compare different neuron
morphology, we use a novel unsupervised graph level representation method to
get the vector representation of each skeleton graphs. We provide detailed
experimental results to demonstrate the effectiveness of our method.
Specifically, our analysis of the Ciona dataset demonstrates that shape could
be used as a significant feature for classifying neuronal types.

##  VI Acknowledgement

The authors would like to thank Kerrianne Ryan for providing the Ciona EM
dataset as well as the great annotations on the dataset. This work was
supported in part by grants from NIH # 5R01NS103774-04 and NSF-SSI award #
1664172.

##  VII Author Contributions

JJ designed the overall neuron morphology analysis pipeline, carried out
experiments on different datasets, and drafted manuscript. MG helped on the
clustering/classification part of the experiments and helped draft the
manuscript. CB provided the Ciona neuron class information and helped draft
the paper. WS coordinated Ciona dataset acquisation, helped analyzed the Ciona
data, and reviewed the manuscript. BSM coordinated the overall design,
development and evaluations of the neuron morphology analysis method, and
helped prepare the manuscript.

## References

  * [1] M. Halavi, K. A. Hamilton, R. Parekh, and G. A. Ascoli, “Digital reconstructions of neuronal morphology: three decades of research trends,” _Frontiers in neuroscience_ , vol. 6, p. 49, 2012.
  * [2] S. Jiang, Z. Pan, Z. Feng, Y. Guan, M. Ren, Z. Ding, S. Chen, H. Gong, Q. Luo, and A. Li, “Skeleton optimization of neuronal morphology based on three-dimensional shape restrictions,” _BMC bioinformatics_ , vol. 21, no. 1, pp. 1–16, 2020.
  * [3] M. Costa, J. D. Manton, A. D. Ostrovsky, S. Prohaska, and G. S. Jefferis, “Nblast: rapid, sensitive comparison of neuronal structure and construction of neuron family databases,” _Neuron_ , vol. 91, no. 2, pp. 293–311, 2016\.
  * [4] M. Abdellah, J. Hernando, S. Eilemann, S. Lapere, N. Antille, H. Markram, and F. Schürmann, “Neuromorphovis: a collaborative framework for analysis and visualization of neuronal morphology skeletons reconstructed from microscopy stacks,” _Bioinformatics_ , vol. 34, no. 13, pp. i574–i582, 2018\.
  * [5] S. Li, T. Quan, C. Xu, Q. Huang, H. Kang, Y. Chen, A. Li, L. Fu, Q. Luo, H. Gong, _et al._ , “Optimization of traced neuron skeleton using lasso-based model,” _Frontiers in neuroanatomy_ , vol. 13, p. 18, 2019.
  * [6] P. K. Saha, Y. Xu, H. Duan, A. Heiner, and G. Liang, “Volumetric topological analysis: a novel approach for trabecular bone classification on the continuum between plates and rods,” _IEEE transactions on medical imaging_ , vol. 29, no. 11, pp. 1821–1838, 2010.
  * [7] O. Panichev and A. Voloshyna, “U-net based convolutional neural network for skeleton extraction,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops_ , 2019, pp. 0–0.
  * [8] T.-C. Lee, R. L. Kashyap, and C.-N. Chu, “Building skeleton models via 3-d medial surface axis thinning algorithms,” _CVGIP: Graphical Models and Image Processing_ , vol. 56, no. 6, pp. 462–478, 1994.
  * [9] S. M. Pizer, J. Hong, J. Vicory, Z. Liu, J. Marron, H.-y. Choi, J. Damon, S. Jung, B. Paniagua, J. Schulz, _et al._ , “Object shape representation via skeletal models (s-reps) and statistical analysis,” in _Riemannian Geometric Statistics in Medical Image Analysis_. Elsevier, 2020, pp. 233–271.
  * [10] C. Lin, C. Li, Y. Liu, N. Chen, Y.-K. Choi, and W. Wang, “Point2skeleton: Learning skeletal representations from point clouds,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2021, pp. 4277–4286.
  * [11] F.-Y. Sun, J. Hoffman, V. Verma, and J. Tang, “Infograph: Unsupervised and semi-supervised graph-level representation learning via mutual information maximization,” in _International Conference on Learning Representations_ , 2019.
  * [12] K. Ryan, Z. Lu, and I. A. Meinertzhagen, “The cns connectome of a tadpole larva of ciona intestinalis (l.) highlights sidedness in the brain of a chordate sibling,” _Elife_ , vol. 5, p. e16962, 2016.
  * [13] L. K. Scheffer, C. S. Xu, M. Januszewski, Z. Lu, S.-y. Takemura, K. J. Hayworth, G. B. Huang, K. Shinomiya, J. Maitlin-Shepard, S. Berg, _et al._ , “A connectome and analysis of the adult drosophila central brain,” _eLife_ , vol. 9, p. e57443, 2020.
  * [14] G. A. Ascoli, D. E. Donohue, and M. Halavi, “Neuromorpho. org: a central resource for neuronal morphologies,” _Journal of Neuroscience_ , vol. 27, no. 35, pp. 9247–9251, 2007.
  * [15] C. R. Qi, L. Yi, H. Su, and L. J. Guibas, “Pointnet++: Deep hierarchical feature learning on point sets in a metric space,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [16] C. R. Qi, H. Su, K. Mo, and L. J. Guibas, “Pointnet: Deep learning on point sets for 3d classification and segmentation,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2017, pp. 652–660.
  * [17] Q.-Y. Zhou, J. Park, and V. Koltun, “Open3d: A modern library for 3d data processing,” _arXiv preprint arXiv:1801.09847_ , 2018.
  * [18] P. V. Tran, “Learning to make predictions on graphs with autoencoders,” in _2018 IEEE 5th international conference on data science and advanced analytics (DSAA)_. IEEE, 2018, pp. 237–245.
  * [19] A. Cardona, S. Saalfeld, J. Schindelin, I. Arganda-Carreras, S. Preibisch, M. Longair, P. Tomancak, V. Hartenstein, and R. J. Douglas, “Trakem2 software for neural circuit reconstruction,” _PloS one_ , vol. 7, no. 6, p. e38011, 2012.
  * [20] T. J. Collins, “Imagej for microscopy,” _Biotechniques_ , vol. 43, no. S1, pp. S25–S30, 2007.
  * [21] C. Yuksel, “Sample elimination for generating poisson disk sample sets,” in _Computer Graphics Forum_ , vol. 34, no. 2. Wiley Online Library, 2015, pp. 25–32.
  * [22] S. Wu, H. Huang, M. Gong, M. Zwicker, and D. Cohen-Or, “Deep points consolidation,” _ACM Transactions on Graphics (ToG)_ , vol. 34, no. 6, pp. 1–13, 2015.
  * [23] M. Kazhdan and H. Hoppe, “Screened poisson surface reconstruction,” _ACM Transactions on Graphics (ToG)_ , vol. 32, no. 3, pp. 1–13, 2013.
  * [24] A. Narayanan, M. Chandramohan, R. Venkatesan, L. Chen, Y. Liu, and S. Jaiswal, “graph2vec: Learning distributed representations of graphs,” _arXiv preprint arXiv:1707.05005_ , 2017.


We differentiate (5.1)1 with respect to $x$ and multiply it by $w_{x}$, and
multiply (5.1)1 by $-z_{xx}$. Finally adding these two equations up, we obtain

$\displaystyle\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}$ $\displaystyle
w_{x}^{2}+\frac{3}{2}\Phi_{3x}w_{x}^{2}+z_{xx}^{2}+z_{x}^{2}+\left\\{\frac{1}{2}\Phi_{3}w_{x}^{2}+\frac{1}{2}ww_{x}^{2}-z_{x}z+(\bar{u}\tilde{u})_{xx}w\right\\}_{x}$
(5.5)
$\displaystyle=-\frac{1}{2}w_{x}^{3}-\Phi_{3xx}ww_{x}+\tilde{u}_{xxx}z_{xx}+(\bar{u}\tilde{u})_{xxx}w.$

Adding (5.4) and (5.5) up, we get

$\displaystyle\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}(w^{2}+w_{x}^{2})+\frac{1}{2}\Phi_{3x}w^{2}+\frac{3}{2}\Phi_{3x}w_{x}^{2}+z_{xx}^{2}+2z_{x}^{2}+z^{2}$
(5.6) $\displaystyle\ \
+\left\\{\frac{1}{2}\Phi_{3}w^{2}+\frac{1}{3}w^{3}+zw+\frac{1}{2}\Phi_{3}w_{x}^{2}+\frac{1}{2}ww_{x}^{2}-2z_{x}z+(\bar{u}\tilde{u})_{xx}w\right\\}_{x}$
$\displaystyle=-\bar{u}\tilde{u}_{x}w-\tilde{u}\bar{u}_{x}w-\tilde{u}_{xxx}z-\frac{1}{2}w_{x}^{3}-\Phi_{3xx}ww_{x}+\tilde{u}_{xxx}z_{xx}+(\bar{u}\tilde{u})_{xxx}w.$

Integrating (5.6) over $\mathbbm{R}_{+}\times(0,t)$, combining it with
$\bar{u}(0)=u_{-}<0$, $\tilde{u}(0,t)=0$ and $w(0,t)=0$, we have

$\displaystyle\|w(t)\|_{1}^{2}+\int_{0}^{t}\left(\|\sqrt{\Phi_{3x}}w(\tau)\|^{2}+\|\sqrt{\Phi_{3x}}w_{x}(\tau)\|^{2}+\|z_{xx}(\tau)\|^{2}+2\|z_{x}(\tau)\|^{2}+\|z(\tau)\|^{2}\right)\,\mathrm{d}\tau-
u_{-}\int_{0}^{t}w_{x}^{2}(0,t)\,\mathrm{d}\tau$ (5.7) $\displaystyle\leq
C\left(\|w_{0}\|_{1}^{2}+\int_{0}^{t}|z_{x}(0,t)||z(0,t)|\,\mathrm{d}\tau+\int_{0}^{t}\int_{\mathbbm{R}_{+}}\left(|\bar{u}|\tilde{u}_{x}|w|+\tilde{u}\bar{u}_{x}|w|+|\tilde{u}_{xxx}|(|z|+|z_{xx}|)\right)\,\mathrm{d}\tau\right.$
$\displaystyle\left.\ \ \ \ \ \ \ \
+\int_{0}^{t}\int_{\mathbbm{R}_{+}}(|w_{x}|^{3}+|\Phi_{3xx}ww_{x}|+|(\bar{u}\tilde{u})_{xxx}w|)\,\mathrm{d}x\mathrm{d}\tau\right).$

Firstly, from (5.1)1, we know that
$z_{x}(0,t)=-u_{-}(w_{x}(0,t)+\tilde{u}_{x}(0,t))$ holds. According to Lemma
2.1 $\mathrm{(ii)}$ with $k=1$, for some small $u_{-}$ corresponding to small
$\delta$ satisfying $Cu_{-}^{2}\leq-\frac{u_{-}}{2}$, we have

$\displaystyle\int_{0}^{t}|z_{x}(0,\tau)||z(0,\tau)|\,\mathrm{d}\tau$
$\displaystyle\leq
C\int_{0}^{t}|z_{x}(0,\tau)|^{2}\,\mathrm{d}\tau+\frac{1}{8}\int_{0}^{t}|z(0,\tau)|^{2}\,\mathrm{d}\tau$
(5.8) $\displaystyle\leq
Cu_{-}^{2}\int_{0}^{t}(|w_{x}(0,\tau)|^{2}+|\tilde{u}_{x}(0,\tau)|^{2})\,\mathrm{d}\tau+\frac{1}{8}\int_{0}^{t}|z(\tau)|_{\infty}^{2}\,\mathrm{d}\tau$
$\displaystyle\leq-\frac{u_{-}}{2}\int_{0}^{t}|w_{x}(0,\tau)|^{2}\,\mathrm{d}\tau+C\delta^{\frac{1}{3}}(1+t)^{-1}+\frac{1}{8}\int_{0}^{t}\|z_{x}(\tau)\|^{2}\,\mathrm{d}\tau+\frac{1}{8}\int_{0}^{t}\|z(\tau)\|^{2}\,\mathrm{d}\tau.$

Secondly, we estimate the right-hand side of (5.7) as follows:

$\int_{\mathbbm{R}_{+}}|\bar{u}|\tilde{u}_{x}|w|\,\mathrm{d}x=\int_{\mathbbm{R}_{+}}(-\bar{u})\tilde{u}_{x}|w|\,\mathrm{d}x=\int_{0}^{u_{+}t}(-\bar{u})\tilde{u}_{x}|w|\,\mathrm{d}x+\int_{u_{+}t}^{\infty}(-\bar{u})\tilde{u}_{x}|w|\,\mathrm{d}x=:I_{1}+I_{2}.$

By virtue of $\bar{u}<0,\tilde{u}>0$, using Lemmas 2.3 and 2.1, we get

$\displaystyle I_{1}$
$\displaystyle\leq|w|_{\infty}|\tilde{u}_{x}|_{\infty}\int_{0}^{u_{+}t}(-\bar{u})\,\mathrm{d}x$
(5.9) $\displaystyle\leq
C\|w_{x}(t)\|^{\frac{1}{2}}\|w(t)\|^{\frac{1}{2}}\delta^{\frac{1}{4}}(1+t)^{-\frac{7}{8}}\int_{0}^{u_{+}t}\frac{\delta}{1+\delta
x}\,\mathrm{d}x$
$\displaystyle\leq\frac{1}{48}\|w_{x}(t)\|^{2}+C\delta^{\frac{1}{3}}\|w(t)\|^{\frac{2}{3}}\left\\{(1+t)^{-\frac{7}{8}}\ln(1+t)\right\\}^{\frac{4}{3}}$
$\displaystyle\leq\frac{1}{48}\|w_{x}(t)\|^{2}+C\delta^{\frac{1}{3}}(1+t)^{-\frac{7}{6}}(\ln(1+t))^{\frac{4}{3}},$

and

$\displaystyle I_{2}$
$\displaystyle\leq|w|_{\infty}\left\\{-[\bar{u}\tilde{u}]_{u_{+}t}^{\infty}+\int_{u_{+}t}^{\infty}\bar{u}_{x}\tilde{u}\,\mathrm{d}x\right\\}$
(5.10) $\displaystyle\leq
Cu_{+}\|w_{x}(t)\|^{\frac{1}{2}}\|w(t)\|^{\frac{1}{2}}\int_{u_{+}t}^{\infty}\bar{u}_{x}\,\mathrm{d}x$
$\displaystyle\leq
C\delta\|w_{x}(t)\|^{\frac{1}{2}}\|w(t)\|^{\frac{1}{2}}\int_{u_{+}t}^{\infty}\frac{\delta^{2}}{(1+\delta
x)^{2}}\,\mathrm{d}x$
$\displaystyle\leq\frac{1}{48}\|w_{x}(t)\|^{2}+C\delta^{\frac{8}{3}}\|w(t)\|^{\frac{2}{3}}(1+t)^{-\frac{4}{3}}$
$\displaystyle\leq\frac{1}{48}\|w_{x}(t)\|^{2}+C\delta^{\frac{1}{3}}(1+t)^{-\frac{4}{3}}.$

In deriving the second inequality of $\eqref{fhbu1}$, we have used from Lemma
2.1 $\mathrm{(iv)}$ and $\mathrm{(v)}$ that

$|\tilde{u}_{x}|_{\infty}\leq
C(\delta(1+t)^{-\frac{1}{2}})^{\frac{1}{4}}((1+t)^{-1})^{\frac{3}{4}}\leq
C\delta^{\frac{1}{4}}(1+t)^{-\frac{7}{8}}.$

And in deriving the last inequality of $\eqref{fhbu1}$ and $\eqref{fhbu2}$, we
have used the fact that $\|w(t)\|$ is bounded. In a similar fashion to (5.9)
and (5.10), we can obtain

$\displaystyle\int_{\mathbbm{R}_{+}}\bar{u}_{x}\tilde{u}|w|\,\mathrm{d}x\leq\frac{1}{24}\|w_{x}\|^{2}+C\delta^{\frac{1}{3}}\left\\{(1+t)^{-\frac{4}{3}}+(1+t)^{-\frac{7}{6}}(\ln(1+t))^{\frac{4}{3}}\right\\}.$
(5.11)

Using Cauchy-Schwarz inequality and Lemma 2.1 $\mathrm{(v)}$ with $k=3,l=0$,
we get

$\int_{\mathbbm{R}_{+}}|\tilde{u}_{xxx}(z+z_{xx})|\,\mathrm{d}x\leq\frac{1}{8}(\|z_{xx}(t)\|^{2}+\|z(t)\|^{2})+C\delta^{\frac{1}{3}}(1+t)^{-\frac{5}{2}},$
(5.12)

and

$\displaystyle\int_{\mathbbm{R}_{+}}|(\bar{u}\tilde{u})_{xxx}w|\,\mathrm{d}x$
$\displaystyle\leq
C\int_{\mathbbm{R}_{+}}(|\bar{u}_{xxx}|\tilde{u}|w|+|\bar{u}_{xx}\tilde{u}_{x}w|+|\bar{u}_{x}\tilde{u}_{xx}w|+|\bar{u}\tilde{u}_{xxx}w|)\,\mathrm{d}x$
(5.13) $\displaystyle\leq
C\left|\frac{\bar{u}_{xxx}}{\bar{u}_{x}}\right|_{\infty}\int_{\mathbbm{R}_{+}}\tilde{u}\bar{u}_{x}|w|\,\mathrm{d}x+C|w|_{\infty}(|\bar{u}_{xx}|_{1}|\tilde{u}_{x}|_{\infty}+|\bar{u}_{x}|_{1}|\tilde{u}_{xx}|_{\infty}+|\bar{u}|_{\infty}|\tilde{u}_{xxx}|_{1})$
$\displaystyle\leq
C\delta|w|_{\infty}\int_{\mathbbm{R}_{+}}\tilde{u}\bar{u}_{x}\,\mathrm{d}x+C\delta|w|_{\infty}(1+t)^{-1}$
$\displaystyle\leq\frac{1}{24}\|w_{x}(t)\|^{2}+C\delta^{\frac{1}{3}}\left\\{(1+t)^{-\frac{4}{3}}+(1+t)^{-\frac{7}{6}}(\ln(1+t))^{\frac{4}{3}}\right\\}.$

From the Lemmas (2.3) and (2.1), we can get estimates of the remaining terms
on the right-hand side of (5.7):

$\displaystyle\int_{\mathbbm{R}_{+}}|w_{x}|^{3}\,\mathrm{d}x\leq|w_{x}|_{\infty}\|w_{x}(t)\|^{2}\leq\varepsilon_{0}\|w_{x}(t)\|^{2},$
(5.14)

and

$\displaystyle\int_{\mathbbm{R}_{+}}|\Phi_{3xx}ww_{x}|\,\mathrm{d}x$
$\displaystyle\leq\int_{\mathbbm{R}_{+}}(|\tilde{u}_{xx}ww_{x}|+|\bar{u}_{xx}ww_{x}|)\,\mathrm{d}x$
(5.15)
$\displaystyle\leq\frac{1}{24}\|w_{x}(t)\|^{2}+C|w|_{\infty}\delta^{\frac{1}{3}}(1+t)^{-\frac{3}{2}}+C\left|\frac{\bar{u}_{xx}^{2}}{\bar{u}_{x}}\right|_{\infty}\|\sqrt{\bar{u}_{x}}w(t)\|^{2}.$

On the other hand, we note that $w_{x}^{2}\leq
3(z_{xx}^{2}+z^{2}+\tilde{u}_{xxx}^{2})$ due to (5.1)2, which yields

$\|w_{x}(t)\|^{2}\leq
3(\|z_{xx}(t)\|^{2}+\|z(t)\|^{2}+\|\tilde{u}_{xxx}(t)\|^{2}).$ (5.16)

In the end, substituting (5.8)-(5.15) into (5.7), and using (5.16), for some
small $\delta$ and $\varepsilon_{0}(<\frac{1}{8})$, we can conclude that

$\displaystyle\|w(t)\|_{1}^{2}+\int_{0}^{t}\left(\|\sqrt{\Phi_{3x}}w(\tau)\|^{2}+\|\sqrt{\Phi_{3x}}w_{x}(\tau)\|^{2}+\|z(\tau)\|_{2}^{2}\right)\,\mathrm{d}\tau-
u_{-}\int_{0}^{t}w_{x}^{2}(0,\tau)\,\mathrm{d}\tau\leq
C(\|w_{0}\|_{1}^{2}+\delta^{\frac{1}{3}}).$

This completes the proof of Lemma 5.4. $\hfill\Box$

From (5.16) and the results of Lemma 5.4, we can easily get the following
estimate.

######  Corollary 5.5.

Under the assumptions of Lemma 5.4, the estimate

$\int_{0}^{t}\|w_{x}(\tau)\|^{2}\,\mathrm{d}\tau\leq
C(\|w_{0}\|_{1}^{2}+\delta^{\frac{1}{3}})$

holds for $t\in[0,T]$.

######  Lemma 5.6.

There are positive constants $\varepsilon_{2}(\leq\varepsilon_{1})$ and
$C=C(\varepsilon_{2})$ such that if $\|w_{0}\|_{2}+\delta\leq\varepsilon_{2}$,
then

$\|w_{xx}(t)\|^{2}+\int_{0}^{t}\left(\|\sqrt{\Phi_{3x}}w_{xx}(\tau)\|^{2}+\|z_{xxx}(\tau)\|^{2}+\|z_{xx}(\tau)\|^{2}+w_{xx}^{2}(0,\tau)\right)\,\mathrm{d}\tau\leq
C(\|w_{0}\|_{2}^{2}+\delta^{\frac{1}{3}}),$ (5.17)

for $t\in[0,T]$.

Proof. We differentiate (5.1)1 twice with respect to $x$ and multiply it by
$w_{xx}$, and differentiate (5.1)2 with respect to $x$ and multiply it by
$-z_{xxx}$. Then, adding these two equations up, we obtain

$\displaystyle\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}w_{xx}^{2}$
$\displaystyle+\frac{5}{2}\Phi_{3x}w_{xx}^{2}+z_{xxx}^{2}+z_{xx}^{2}+\left\\{\frac{1}{2}\Phi_{3}w_{xx}^{2}+\frac{1}{2}ww_{xx}^{2}-z_{xx}z_{x}+(\bar{u}\tilde{u})_{xxx}w_{x}\right\\}_{x}$
(5.18)
$\displaystyle=-\Phi_{3xxx}ww_{xx}-3\Phi_{3xx}w_{x}w_{xx}-\frac{5}{2}w_{x}w_{xx}^{2}+\partial_{x}^{4}\tilde{u}z_{xxx}+\partial_{x}^{4}(\bar{u}\tilde{u})w_{x}.$

Integrating (5.18) over $\mathbbm{R}_{+}\times(0,t)$, combining
$\Phi_{3}(0,t)=u_{-}<0$ and $w(0,t)=0$, we have

$\displaystyle\|w_{xx}(t)\|^{2}+\int_{0}^{t}\left(\|\sqrt{\Phi_{3x}}w_{xx}(\tau)\|^{2}+\|z_{xxx}(\tau)\|^{2}+\|z_{xx}(\tau)\|^{2}\right)\,\mathrm{d}\tau-
u_{-}\int_{0}^{t}w_{xx}^{2}(0,\tau)\,\mathrm{d}\tau$ (5.19) $\displaystyle\leq
C\left(\|w_{0}\|_{2}^{2}+\int_{0}^{t}|z_{xx}(0,\tau)||z_{x}(0,\tau)|\,\mathrm{d}\tau+\int_{0}^{t}|(\bar{u}\tilde{u})_{xxx}(0,\tau)||w_{x}(0,\tau)|\,\mathrm{d}\tau\right.$
$\displaystyle\left.\ \ \ \
+\int_{0}^{t}\int_{\mathbbm{R}_{+}}(|\Phi_{3xxx}ww_{xx}|+|\Phi_{3xx}w_{x}w_{xx}|+|w_{x}w_{xx}^{2}|+|\partial_{x}^{4}\tilde{u}z_{xxx}|+|\partial_{x}^{4}(\bar{u}\tilde{u})w_{x}|)\,\mathrm{d}x\mathrm{d}\tau\right).$

Firstly, we estimate the second and the third terms on the right-hand side of
(5.19). Using the Cauchy-Schwarz inequality, we have

$\displaystyle C\int_{0}^{t}|z_{xx}(0,\tau)||z_{x}(0,\tau)|\,\mathrm{d}\tau$
$\displaystyle\leq
C\int_{0}^{t}|z_{xx}(\tau)|_{\infty}|z_{x}(\tau)|_{\infty}\,\mathrm{d}\tau$
(5.20)
$\displaystyle\leq\frac{1}{4}\int_{0}^{t}\|z_{xxx}(\tau)\|^{2}\,\mathrm{d}\tau+C\int_{0}^{t}(\|z_{xx}(\tau)\|^{2}+\|z_{x}(\tau)\|^{2})\,\mathrm{d}\tau$
$\displaystyle\leq\frac{1}{4}\int_{0}^{t}\|z_{xxx}(\tau)\|^{2}\,\mathrm{d}\tau+C(\|w_{0}\|_{1}^{2}+\delta^{\frac{1}{3}}).$

For
$(\bar{u}\tilde{u})_{xxx}=\bar{u}_{xxx}\tilde{u}+3\bar{u}_{xx}\tilde{u}_{x}+3\bar{u}_{x}\tilde{u}_{xx}+\bar{u}\tilde{u}_{xxx}$
and $\tilde{u}(0,t)=0$, combining Lemma 5.4 and Lemma 2.1 $\mathrm{(ii)}$, we
get

$\displaystyle
C\int_{0}^{t}|(\bar{u}\tilde{u})_{xxx}(0,\tau)||w_{x}(0,\tau)|\,\mathrm{d}\tau$
$\displaystyle\leq
C\delta\int_{0}^{t}(|\tilde{u}_{x}(0,\tau)|+|\tilde{u}_{xx}(0,\tau)|+|\tilde{u}_{xxx}(0,\tau)|)|w_{x}(0,\tau)|\,\mathrm{d}\tau$
$\displaystyle\leq\int_{0}^{t}w_{x}^{2}(0,\tau)\,\mathrm{d}\tau+C\delta^{\frac{1}{3}}(1+t)^{-1}.$

Next, we estimate the last five terms on the right-hand side of (5.19). Using
Cauchy-Schwarz inequality, we have

$\displaystyle\int_{\mathbbm{R}_{+}}|\Phi_{3xxx}ww_{xx}|\,\mathrm{d}x$
$\displaystyle\leq\int_{\mathbbm{R}_{+}}(|\bar{u}_{xxx}ww_{xx}|+|\tilde{u}_{xxx}ww_{xx}|)\,\mathrm{d}x$
$\displaystyle\leq\frac{1}{36}\|w_{xx}(t)\|^{2}+C\left|\frac{\bar{u}_{xxx}^{2}}{\bar{u}_{x}}\right|_{\infty}\|\sqrt{\bar{u}_{x}}w(t)\|^{2}+C\delta^{\frac{1}{3}}|w(t)|_{\infty}(1+t)^{-\frac{5}{2}},$
$\displaystyle\int_{\mathbbm{R}_{+}}|\Phi_{3xx}w_{x}w_{xx}|\,\mathrm{d}x\leq\frac{1}{36}\|w_{xx}(t)\|^{2}+|\Phi_{3xx}(t)|_{\infty}\|w_{x}(t)\|^{2}\leq\frac{1}{36}\|w_{xx}(t)\|^{2}+C\delta\|w_{x}(t)\|^{2},$
$\displaystyle\int_{\mathbbm{R}_{+}}|\partial_{x}^{4}\tilde{u}z_{xxx}|\,\mathrm{d}x\leq\frac{1}{4}\|z_{xxx}(t)\|^{2}+C\delta^{\frac{1}{3}}(1+t)^{-\frac{7}{2}},$

and by utilizing the _a priori_ assumption (5.3),

$\displaystyle\int_{\mathbbm{R}_{+}}|w_{x}w_{xx}^{2}|\,\mathrm{d}x\leq|w_{x}(t)|_{\infty}\|w_{xx}(t)\|^{2}\leq\varepsilon_{0}\|w_{xx}(t)\|^{2}.$

By the same method as (5.9) and (5.10), we can get the last one estimate:

$\displaystyle\int_{\mathbbm{R}_{+}}|\partial_{x}^{4}(\bar{u}\tilde{u})w_{x}|\,\mathrm{d}x$
$\displaystyle\leq
C\int_{\mathbbm{R}_{+}}(|\partial_{x}^{4}\bar{u}\tilde{u}|+|\bar{u}_{xxx}\tilde{u}_{x}|+|\bar{u}_{xx}\tilde{u}_{xx}|+|\bar{u}_{x}\tilde{u}_{xxx}|+|\bar{u}\partial_{x}^{4}\tilde{u}|)|w_{x}|\,\mathrm{d}x$
(5.21) $\displaystyle\leq
C\left|\frac{\partial_{x}^{4}\bar{u}}{\bar{u}_{x}}\right|_{\infty}|w_{x}|_{\infty}\int_{\mathbbm{R}_{+}}\bar{u}_{x}\tilde{u}\,\mathrm{d}x+C\delta|w_{x}|_{\infty}(|\tilde{u}_{x}|_{\infty}+|\tilde{u}_{xx}|_{\infty}+|\tilde{u}_{xxx}|_{\infty}+|\partial_{x}^{4}\tilde{u}|_{\infty})$
$\displaystyle\leq
C\delta|w_{x}|_{\infty}\left(\int_{0}^{u_{+}t}\bar{u}_{x}\tilde{u}\,\mathrm{d}x+\int_{u_{+}t}^{\infty}\bar{u}_{x}\tilde{u}\,\mathrm{d}x\right)+C\delta|w_{x}|_{\infty}(1+t)^{-1}$
$\displaystyle\leq\frac{1}{36}\|w_{xx}(t)\|^{2}+C\delta^{\frac{1}{3}}\|w_{x}(t)\|^{\frac{2}{3}}\left\\{(1+t)^{-\frac{4}{3}}+(1+t)^{-\frac{7}{6}}(\ln(1+t))^{\frac{4}{3}}\right\\}.$

In the end, to estimate $\int_{0}^{t}\|w_{xx}(\tau)\|^{2}\,\mathrm{d}\tau$,
for (5.1)2, we note that

$\|w_{xx}(t)\|^{2}\leq
3(\|z_{xxx}(t)\|^{2}+\|z_{x}(t)\|^{2}+\|\partial_{x}^{4}\tilde{u}(t)\|^{2}).$
(5.22)

Substituting (5.20)-(5.21) into (5.19), combining (5.22) and Lemma 5.4, we
have

$\displaystyle\|w_{xx}(t)\|^{2}+\int_{0}^{t}\left(\|\sqrt{\Phi_{3x}}w_{xx}(\tau)\|^{2}+\|z_{xxx}(\tau)\|^{2}+\|z_{xx}(\tau)\|^{2}+w_{xx}^{2}(0,\tau)\right)\,\mathrm{d}\tau$
$\displaystyle\leq$
$\displaystyle\left(\frac{3}{4}+\varepsilon_{0}\right)\|z_{xxx}(t)\|^{2}+C(\|w_{0}\|_{2}^{2}+\delta^{\frac{1}{3}}).$

That is, for $\varepsilon_{0}<\frac{1}{4}$, we complete the proof of (5.17).
$\hfill\Box$

The equation (5.22) and Lemma 5.6 yields Corollary 5.7.

######  Corollary 5.7.

Under the assumptions of Lemma 5.6,

$\int_{0}^{t}\|w_{xx}(\tau)\|^{2}\,\mathrm{d}\tau\leq
C(\|w_{0}\|_{2}^{2}+\delta^{\frac{1}{3}}),\quad\forall t\in[0,T].$

######  Lemma 5.8.

Under the assumptions of Lemma 5.6,

$\|z(t)\|_{3}^{2}\leq C(\|w_{0}\|_{2}^{2}+\delta^{\frac{1}{3}}),\quad\forall
t\in[0,T].$ (5.23)

Proof. Firstly, from (5.1)1, we have the equation at the boundary $x=0$,

$z_{x}(0,t)=-u_{-}(w_{x}(0,t)+\tilde{u}_{x}(0,t)),$

which plays an essential role in estimating boundary terms. From (5.1)2, it
holds

$z_{xx}^{2}+z^{2}+2z_{x}^{2}=w_{x}^{2}+\tilde{u}_{xxx}^{2}+2(z_{x}z)_{x}+2w_{x}\tilde{u}_{xxx}.$
(5.24)

Integrating (5.24) over $\mathbbm{R}_{+}$, by Cauchy-Schwarz inequality, we
get

$\displaystyle\|z_{xx}(t)\|^{2}+2\|z_{x}(t)\|^{2}+\|z(t)\|^{2}$
$\displaystyle\leq
2\|w_{x}(t)\|^{2}+2\|\tilde{u}_{xxx}(t)\|^{2}+2|z_{x}(0,t)||z(0,t)|$ (5.25)
$\displaystyle\leq
C(\|w_{0}\|_{2}^{2}+\delta^{\frac{1}{3}})+C|z_{x}(0,t)|^{2}+\frac{1}{2}|z(t)|_{\infty}^{2}$
$\displaystyle\leq
C(\|w_{0}\|_{2}^{2}+\delta^{\frac{1}{3}})+C(|w_{x}(t)|_{\infty}^{2}+|\tilde{u}_{x}(t)|_{\infty}^{2})+\frac{1}{2}\|z_{x}(t)\|^{2}+\frac{1}{2}\|z(t)\|^{2}$
$\displaystyle\leq
C(\|w_{0}\|_{2}^{2}+\delta^{\frac{1}{3}})+\frac{1}{2}\|z_{x}(t)\|^{2}+\frac{1}{2}\|z(t)\|^{2}.$

Differentiating (5.1)2 with respect to $x$ and integrating the resulting
equation over $\mathbbm{R}_{+}$, combining (5.25), Lemma 2.1 and Lemma 5.6, we
get

$\displaystyle\|z_{xxx}(t)\|^{2}\leq
C(\|z_{x}(t)\|^{2}+\|w_{xx}(t)\|^{2}+\|\partial_{x}^{4}\tilde{u}(t)\|^{2})\leq
C(\|w_{0}\|_{2}^{2}+\delta^{\frac{1}{3}}).$ (5.26)

The combination of (5.25) and (5.26) completes the proof of Lemma 5.8.
$\hfill\Box$

###  5.3 Asymptotic Behavior toward the Superposition of Nonlinear Waves

Once the _a priori_ estimates is established, by combining the local
existence, the global existence of unique solution of (5.1) and its asymptotic
behavior are easily obtained. That is, the global in time solution

$\begin{cases}w\in C^{0}([0,\infty);H^{2}),\quad w_{x}\in
L^{2}(0,\infty;H^{1}),\\\\[2.84526pt] z\in C^{0}([0,\infty);H^{3})\cap
L^{2}(0,\infty;H^{3}).\end{cases}$

Then, the _a priori_ estimates again assert that

$\begin{cases}\sup\limits_{t\geq
0}(\|w(t)\|_{2}^{2}+\|z(t)\|_{3}^{2})<\infty,\\\\[2.84526pt]
\int_{0}^{t}(\|w_{x}\|_{1}^{2}+\|z\|_{3}^{2})\,\mathrm{d}\tau<\infty.\end{cases}$
(5.27)

To complete the proof of Theorem 5.1, we need to show that

$\int_{0}^{\infty}\left|\frac{\mathrm{d}}{\mathrm{d}t}\|w_{x}\|^{2}\right|\,\mathrm{d}t<\infty,$
(5.28)

it follows

$\|w_{x}\|\rightarrow 0,\ \ (t\rightarrow 0).$ (5.29)

Using the Sobolev inequality, we can obtain the desired asymptotic behavior in
Theorem 5.1

$\begin{cases}\sup\limits_{x\in\mathbbm{R}_{+}}|w(x,t)|\leq\sqrt{2}\|w\|^{\frac{1}{2}}\|w_{x}\|^{\frac{1}{2}}\rightarrow
0,\ \ (t\rightarrow 0),\\\\[11.38109pt]
\sup\limits_{x\in\mathbbm{R}_{+}}|w_{x}(x,t)|\leq\sqrt{2}\|w_{x}\|^{\frac{1}{2}}\|w_{xx}\|^{\frac{1}{2}}\rightarrow
0,\ \ (t\rightarrow 0).\end{cases}$ (5.30)

The combination of (5.29) and (5.30) can completes the proof of Theorem 5.1.

The proof of (5.28) can be easily obtained. In fact, using the similar
estimatea as $\eqref{eq-wx3}$, $\eqref{fhbwx3}$, $\eqref{eq-uxxxzxx}$ and
$\eqref{eq-uuw}$, combining (5.27), we can get from $\eqref{fhbwx1}$ that

$\displaystyle\int_{0}^{\infty}\left|\frac{\mathrm{d}}{\mathrm{d}t}\|w_{x}\|^{2}\right|\,\mathrm{d}t$
$\displaystyle\leq
C\int_{0}^{\infty}|z_{x}(0,t)z(0,t)|\,\mathrm{d}t+C\int_{0}^{\infty}\int_{\mathbbm{R}_{+}}(|w_{x}^{3}|+|\Phi_{3xx}ww_{x}|+|\tilde{u}_{xxx}z_{xx}|+|(\bar{u}\tilde{u})_{xxx}w|)\,\mathrm{d}x\mathrm{d}t$
$\displaystyle<\infty.$

Thus, we get the asymptotic behavior of $w$ and we finish the proof of (5.2)1
in Theorem 5.1. In the end, using the Lemma 2.12, and setting
$g=-u_{-}(w_{x}(0,t)+\tilde{u}_{x}(0,t))$,
$f(x)=-w_{x}(x,t)-\tilde{u}_{xxx}(x,t)$ for any fixed $t\in[0,\infty)$, we can
obtain the asymptotic behavior of $z$, which completes the proof of Theorem
5.1.

Acknowledgements: The research was supported by Guangdong Basic and Applied
Basic Research Foundation $\\#$2020B1515310015, $\\#$2021A1515010367, the
National Natural Science Foundation of China $\\#$11771150, $\\#$11831003.

## References

  * [1] M. Di Francesco, Initial value problem and relaxation limits of the Hamer model for radiating gases in several space variables, Nonlinear Differ. Equ. Appl., 13(2007), 531-562.
  * [2] M. Di Francesco, C. Lattanzio, Optimal $L^{1}$ rate of decay to diffusion waves for the Hamer model of radiating gases, Appl. Math. Lett., 19(2006), 1046-1052.
  * [3] R.J. Duan, K. Fellner, C.J. Zhu, Energy method for multi-dimensional balance laws with non-local dissipation, J. Math. Pures Appl., 93(2010), 572-598.
  * [4] R.J. Duan, L.Z. Ruan, C.J. Zhu, Optimal decay rates to conservation laws with diffusion-type terms of regularity-gain and regularity-loss, Math. Models Methods Appl. Sci., 22(2012), 1250012, 1-39.
  * [5] W.L. Gao, L.Z. Ruan, C.J. Zhu, Decay rates to the planar rarefaction waves for a model system of the radiating gas in $n$ dimensions, J. Differential Equations, 244(2008), 2614-2640.
  * [6] W.L. Gao, C.J. Zhu, Asymptotic decay toward the planar rarefaction waves for a model system of the radiating gas in two dimensions, Math. Models Methods Appl. Sci., 18(2008), 511-541.
  * [7] K. Hamer, Nonlinear effects on the propagation of sound waves in a radiating gas, Quart. J. Mech. Appl. Math., 24(1971), 155-168.
  * [8] I. Hashimoto, A. Matsumura, Large time behavior of solutions to an initial boundary value problem on the half space for scalar viscous conservation law, Methods Appl. Anal., 14(2007), 45-59.
  * [9] Y. Hattori, K. Nishihara, A note on the stability of rarefaction wave of the Burgers equation, Jpn. J. Ind. Appl. Math., 8(1991), 85-96.
  * [10] M.A. Heaslet, B.S. Baldwin, Predictions of the structure of radiation resisted shock waves, Phys. Fluids, 6(1963),781-791.
  * [11] L. Hsiao, T.-P. Liu, Convergence to nonlinear diffusion waves for solutions of a system of hyperbolic conservation laws with damping, Comm. Math. Phys., 143(1992), 599-605.
  * [12] F.M. Huang, R.H. Pan, Convergence rate for compressible Euler equations with damping and vacuum, Arch. Ration. Mech. Anal., 166(2003), 359-376.
  * [13] K. Ito, Asymptotic decay toward the planar rarefaction waves for viscous conservation laws in several dimentions, Math. Mod. Meth. Appl. Sci., 6(1996), 315-338.
  * [14] S. Kawashima, S. Nishibata, Shock waves for a model system of a radiating gas, SIAM J. Math. Anal., 30(1999), 95-117.
  * [15] S. Kawashima, Y. Tanaka, Stability of rarefaction waves for a model system of a radiating gas, Kyushu J. Math., 58(2004), 211-250.
  * [16] C. Lattanzio, P. Marcati, Global well-posedness and relaxation limits of a model for radiating gas, J. Differential Equations, 190(2003), 439-465.
  * [17] H. Liu, E. Tadmor, Critical thresholds in a conservation model for nonlinear conservation laws, SIAM J. Math. Anal., 33(2001), 930-945.
  * [18] T.-P. Liu, A. Matsumura, K. Nishihara, Behaviors of solutions for the Burgers equation with boundary corresponding to rarefaction waves, SIAM J. Math. Anal., 29(1998), 293-308.
  * [19] T.-P. Liu, K. Nishihara, Asymptotic behavior for scalar viscous conservation laws with boundary effect, J. Differential Equations, 133(1997), 296-320.
  * [20] T.-P. Liu, S.-H. Yu, Propagation of a stationary shock layer in the presence of a boundary, Arch. Ration. Mech. Anal., 139(1997), 57-82.
  * [21] Y.Q. Liu, S. Kawashima, Asymptotic behavior of solutions to a model system of a radiating gas, Comm. Pure Appl. Anal., 10(2011), 209-223.
  * [22] R.B. Lowrie, J.D. Edwards, Radiative shock solutions with grey nonequilibrium diffusion, Shock Waves, 18(2008), 129-143.
  * [23] R.B. Lowrie, R.M. Rauenzahn, Radiative shock solutions in the equilibrium diffusion limit, Shock Waves, 16(2007), 445-453.
  * [24] T. Luo, Asymptotic stability of planar rarefaction waves for relaxation approximation of conservation laws in several dimensions, J. Differential Equations, 133(1997), 255-279.
  * [25] A. Matsumura, M. Mei, Convergence to travelling fronts of solutions of the p-system with viscosity in the presence of a boundary, Arch. Ration. Mech. Anal., 146(1999), 1-22.
  * [26] A. Matsumura, K. Nishihara, Global stability of the rarefaction waves of a one-dimensional model system for compressible viscous gas, Comm. Math. Phys., 144(1992), 325-335.
  * [27] T. Nakamura, Asymptotic decay toward the rarefaction waves of solutions for viscous conservation laws in a one dimensional half space, SIAM J. Math. Anal., 34(2003), 1308-1317.
  * [28] T. Nguyen, Multi-dimensional stability of planar Lax shocks in hyperbolic-elliptic coupled systems, J. Differential Equations, 252(2012), 382-411.
  * [29] M. Nishikawa, K. Nishihara, Asymptotics toward the planar rarefaction wave for viscous conservation law in two space dimensions, Trans. Amer. Math. Soc., 352(2000), 1203-1215.
  * [30] M. Ohnawa, $L^{\infty}$-stability of continuous shock waves in a radiating gas model, SIAM J. Math. Anal., 46(2014), 2136-2159.
  * [31] L.Z. Ruan, J. Zhang, Asymptotic stability of rarefaction wave for hyperbolic-elliptic coupled system in radiating gas, Acta Math. Sci. Ser. B Engl. Ed., 27(2007) 347-360.
  * [32] L.Z. Ruan, C.J. Zhu, Asymptotic decay toward rarefaction wave for a hyperbolic-elliptic coupled system on half space, J. Partial Differential Equations, 21(2008), 173-192.
  * [33] L.Z. Ruan, C.J. Zhu, Asymptotic behavior of solutions to a hyperbolic-elliptic coupled system in multi-dimensional radiating gas, J. Differential Equations, 249(2010), 2076-2110.
  * [34] S. Schochet, E. Tadmor, The regularized Chapman-Enskog expansion for scalar conservation laws, Arch. Ration. Mech. Anal., 119(1992), 95-107.
  * [35] D. Serre, $L^{1}$-stability of constants in a model for radiating gases, Comm. Math. Sci., 1(2003), 197-205.
  * [36] Y. Tanaka, Asymptotic behavior of solutions to the one-dimensional model system for a radiating gas, Master’s Thesis, Kyushu University, 1995 (in Japanese).
  * [37] W. G. Vincenti, C.H. Kruger, Introduction to Physical Gas Dynamics, Wiley, New York, 1965.
  * [38] W.J. Wang, W.K. Wang, Z.G. Wu, Decay of a model system of radiating gas, Math. Methods Appl. Sci., 37(2014), 2331-2340.
  * [39] Z.P. Xin, Asymptotic stability of rarefaction waves for $2\times 2$ viscous hyperbolic conservation laws-the two modes case, J. Differential Equations, 78(1989), 191-219.
  * [40] T. Yang, H.J. Zhao, BV estimates on Lax-Friedrichs’ scheme for a model of radiating gas, Appl. Anal., 83(2004), 533-539.
  * [41] T. Yang, H.J. Zhao, C.J. Zhu, Asymptotic behavior of solution to a hyperbolic system with relaxition and boundary effect, J. Differential Equations, 163(2000), 348-380.
  * [42] Ya.B. Zel’dovich, Shock waves of large amplitude in air, Soviet Phys. JETP, 5(1957), 919-927.
  * [43] Ya.B. Zel’dovich, Yu.P. Raizer, Physics of Shock Waves and High-temperature Hydrodynamic Phenomena, Academic Press, 1966.
  * [44] H.J. Zhao, Nonlinear stability of strong planar rarefaction waves for the relaxation approximation of conservation laws in several space dimensions, J. Differential Equations, 163(2000), 198-223.
  * [45] C.J. Zhu, Asymptotic behavior of solutions for p-system with relaxation, J. Differential Equations, 180(2002), 273-306.