#  _MAROAM_ : Map-based Radar SLAM through Two-step Feature Selection

Dequan Wang, Yifan Duan, Xiaoran Fan, Chengzhen Meng, Jianmin Ji and Yanyong
Zhang School of Computer Science and Technology, University of Science and
Technology of China, Hefei, 230026, China {wdq15588, dyf0202,
<EMAIL_ADDRESS><EMAIL_ADDRESS>{jianmin,
<EMAIL_ADDRESS>
###### Abstract

In this letter, we propose _MAROAM_ , a millimeter wave radar-based SLAM
framework, which employs a two-step feature selection process to build the
global consistent map. Specifically, we first extract feature points from raw
data based on their local _geometric_ properties to filter out those points
that violate the principle of millimeter-wave radar imaging. Then, we further
employ another round of _probabilistic_ feature selection by examining how
often and how recent the feature point has been detected in the proceeding
frames. With such a two-step feature selection, we establish a global
consistent map for accurate and robust pose estimation as well as other
downstream tasks. At last, we perform loop closure and graph optimization in
the back-end, further reducing the accumulated drift error.

We evaluate the performance of _MAROAM_ on the three datasets: the Oxford
Radar RobotCar Dataset, the MulRan Dataset and the Boreas Dataset. We consider
a variety of experimental settings with different scenery, weather, and road
conditions. The experimental results show that the accuracy of _MAROAM_ is
7.95%, 37.0% and 8.9% higher than the currently best-performing algorithms on
these three datasets, respectively. The ablation results also show that our
map-based odometry performs 28.6% better than the commonly used scan-to-frames
method. Finally, as devoted contributors to the open-source community, we will
open source the algorithm after the paper is accepted.

##  I INTRODUCTION

In recent years, millimeter wave radar (radar, in short) has received an
increasing amount of attention in autonomous systems. Compared to vision and
LiDAR, radar offers longer detection ranges and more robust performance under
adverse lighting and weather conditions, lending it towards a viable sensing
option for applications like self driving and robotics[1]. However, in order
for radar to grow from a supplementary sensing approach to one of the primary
ones, their sensing capability, especially when used alone, needs to be
carefully investigated and strengthened. Towards this goal, we study how to
utilize the radar alone for the Simultaneous Localization and Mapping (SLAM)
task in this work.

The main weakness for millimeter wave radar images lies in the fact that they
suffer from the multipath effect, sidelobes, clutter, etc.[1], which may lead
to false alarms. In order to overcome the effect of false alarms, the recent
radar SLAM systems mainly focus on extracting robust features from radar
images. For example, Cen et al. proposed to extract features using traditional
signal processing techniques[2, 3], while Barnes et al. proposed end-to-end
learning-based feature extraction methods[4, 5]. However, even with these
efforts, feature extraction for radar systems remains an active research
topic, which aims at minimizing the cumulative errors in the feature matching
task.

Figure 1: An example radar point cloud map with global consistency, which is
built by _MAROAM_ with the Oxford Radar RobotCar Dataset. For clarity, we also
zoom in several locations and show the detailed views. In this example, we use
the ”17-13-26” sequence with a total length of $10km$.

Next, let us take a close look at the radar odometry. When the scan-to-map
registration strategy is engaged, the number of features grows rapidly with
the map, a large portion of which are of low quality and should not be
selected. The quality of the map seriously affects the performance of map
based-SLAM[6]. RadarSLAM[7] uses descriptor-based pairwise matching of
features. However, since image noise and different radar cross-sectional (RCS)
generated by different viewing angles[8] will directly affect the
representation of descriptors, so it has high requirements for the extraction
process of features and descriptors.

Since pairwise matching of descriptors is not required, the Iterative Closest
Point (ICP)[9] methods achieve state-of-the-art accuracy in motion estimation,
such as CFEAR RadarOdometry[10] and NDT[11], but the submaps they use are only
the superposition of multiple frames, and there is no mapping process. The
disadvantage of ICP-based point cloud registration algorithms is that they are
sensitive to initial values and tend to converge to local optimum[12]. Since
radar point cloud has more noise, the ICP-based algorithm is more likely to
fall into a local optimum, resulting in a larger motion estimation error. We
have demonstrated this experimentally at IV-C. Considering the importance of
map quality, our paper looks at how feature selection can be used to optimize
the process of map building and scan-to-map matching in ICP-based radar SLAM,
resulting in more accurate maps and more robust state estimates.

In order to address these limitations, we propose a Map-based Radar Odometry
and Mapping system, _MAROAM_ , which is an accurate and robust map-based radar
SLAM framework based on LOAM[13]. In order to improve map quality, we use a
two-step geometry-probability feature selection strategy in the registration
and map update stage. In specificity, we first extract surface features by
calculating the Local Linearity and the Local Aggregation with a _geometry-
based_ filter. We then involve a point-to-line ICP method, using the surface
features to solve the relative pose transformation. Subsequently, we use a
_probability-based_ filter to dynamically filter out the features that are not
frequently selected in the ICP matching process. Using the more frequent
features, our algorithm realizes the scan-to-map matching strategy of ICP-
based radar SLAM for the first time. Additionally, We use the Scan Context[14]
for loop closure and modify it to fit the radar data format.

We evaluate the performance of _MAROAM_ on the three datasets: the Oxford
Radar RobotCar Dataset[7], the MulRan Dataset[15] and the Boreas Dataset[16].
We consider a variety of experimental settings with different scenery,
weather, and road conditions. The experimental results show that the accuracy
of _MAROAM_ is 7.95%, 37.0% and 8.9% higher than the currently best-performing
algorithms on these three datasets, respectively. The ablation results also
show that our map-based odometry performs 28.6% better than the commonly used
scan-to-frames method.

In summary, our main contributions are as follows:

  * •

We propose _MAROAM_ , an map-based radar SLAM framework based on LOAM. Using a
geometry-probability two-step feature selection, our algorithm realizes the
scan-to-map matching strategy of ICP-based radar SALM for the first time.

  * •

We evaluate our _MAROAM_ on three datasets which contain a variety of
scenarios, weathers, and road conditions. Experiments show that our algorithm
outperforms SOTA algorithm on all three datasets.

##  II Related Works

###  II-A Radar SLAM

We divide the existing radar odometry/SLAM algorithms into four categories
according to the method of estimating the relative pose:

####  II-A1 Direct methods

Checchin et al. [17] first time use the Fourier-Mellin Transform (FMT) to
register radar images in a sequence for motion estimation. PhaRao[18] apply
FMT to Cartesian and log polar radar images to estimate rotation and
translation to decouple rotation and translation.

####  II-A2 ICP-based methods

Adolfsson et al.[10] process the landmarks extracted by CFAR through
k-strongest filtering in Cartesian coordinates, and optimizes the point-to-
line optimize metric to estimate the relative pose. Kung et al.[19] propose a
RO with probabilistic submap building, and an NDT-based radar scan matching.
Both of them use scan-to-submap matching, but the submaps they use are only
the superposition of multiple frames, and there is no mapping process.

####  II-A3 Descriptor-based methods

Considering the unnecessary influence of radar echo, Cen et al.[2] propose an
algorithm to extract landmarks, and performe scan matching by greedily adding
features correspondence based on unary descriptor and pairwise compatibility
score. After that, another feature extraction algorithm that only use one
parameter is proposed by Cen et al.[3], and graph matching is used for scan
matching. Hone et al.[20, 21] propose a full radar-based SLAM pipeline,
RadarSLAM, composed of pose tracking, local mapping, loop closure detection
and pose graph optimization. RadarSLAM uses visual features for scan matching
and M2DP[22] descriptor based on point cloud for loop detection.

####  II-A4 Learning-based methods

Barnes et al.[5] use deep neural network to learn an embedding space that is
basically free of artifacts and interference, which is used to perform
effective correlation matching between continuous radar scans, and achieve
high accuracy without considering spatial cross-validation(CSV). A self
supervised learning framework is proposed by Barnes et al.[4] to detect the
robust key features of range estimation and metric positioning in radar.
Burnett et al.[23] uses unsupervised method to extract features, which
increases generalization and performs well on both Oxford and Boreas.

###  II-B LiDAR SLAM

LOAM[13] is the pioneering work of 3D LiDAR SLAM, which proposes a basic
framework. LOAM extracts edge features and plane features from the original
point cloud, and designs related loss functions for each type of feature. The
matching consists of a fast frame-to-frame match and a slow frame-to-graph
match. F-LOAM[24] follows LOAM and abandons the scan-to-scan match and
replaces it by only scan-to-map with high frequency. Duan et al.[25] propose a
feature filter, _pFilter_ , by properly measuring each feature point’s p-Index
and only keeping those with high index values, and improved both the
efficiency and accuracy of the registration process. Our work is inspired by
the aforementioned LiDAR SLAM.

Figure 2: Pipeline of the _MAROAM_. The raw point cloud represents the point
cloud directly extracted by the feature extractor, point cloud surf represents
the surface features filtered by _geometry-based_ filter, raw local map is the
union of the local map and the current point cloud surf, local map is
generated from raw local map filtered by _probability-based_ filter.

##  III Method

(a) raw point cloud

(b) point cloud surf

(c) Local map before
probability-based filter

(d) Local map after
probability-based filter

Figure 3: The effect of _geometry-based_ filter and _probability-based_ filter
for a radar scan and a local feature map. The raw point cloud (a) contain a
large number of wrong feature points such as ghosts and speckles; (b) shows
the point cloud after _geometry-based_ filter; (c) shows the local surface map
by stacking all surface feature points; (d) shows local surface map after
_probability-based_ filter. Each feature point in the maps is colored by its
hit probability value: the red color indicates feature points with low hit
probability values, i.e., transient feature points, while the green color
indicates feature points with high hit probability values, i.e., frequent
feature points. The resulting maps in (d) contain fewer feature points that
should be included for registration.

###  III-A System Overview

An overview of _MAROAM_ is shown in Fig. 2. The system receives input from a
2D radar image and outputs a 3-DOF pose estimation and a global consistent
map. The SLAM system consists of two modules: tracking and loop closure. The
tracking module is used for feature extraction, feature filtering and real-
time motion estimation, and the loop closure module is used for loop detection
and graph optimization.

###  III-B Geometry-based Feature Filter

Cen2019[3] is used as feature detector in _MAROAM_. Feature detector firstly
identifies continuous regions with high intensity and low gradient for each
azimuth of a radar image. Then, the midpoints of each continuous region are
extracted as features. The point cloud extracted from the kth scan is denoted
as $\mathcal{P}_{k}$.

The _geometry-based_ filter is used to extract surface features from raw point
cloud, marked as $\mathcal{S}_{k}$. Surface features are defined as points
with high aggregation and high local linearity. Firstly, we use the principal
components analysis (PCA) to calculate the the local linearity. For each
feature $\textbf{p}_{i}\in\mathcal{P}_{k}$, we find the nearest m points to
$\textbf{p}_{i}$ from the point cloud $\mathcal{P}_{k}$ to form a
corresponding set, marked $\mathcal{N}_{i}$. We store the point cloud
$\mathcal{P}_{k}$ in a KD-tree to accelerate nearest neighbor search. Let
$\bar{\textbf{p}}_{i}$ be the centroid of $\mathcal{N}_{i}$ , and
$\mathbf{C}_{i}$ be the covraiance matrix of $\mathcal{N}_{i}$ which is
calculated as

$\mathbf{C}_{i}=\frac{1}{|\mathcal{N}_{i}|}\sum_{\textbf{p}_{j}\in\mathcal{N}_{i}}(\textbf{p}_{j}-\bar{\textbf{p}}_{i})((\textbf{p}_{j}-\bar{\textbf{p}}_{i})^{T},$
(1)

Then We use PCA to calculate the eigenvalues of matrix $\mathbf{C}_{i}$. The
largest two eigenvalues are recorded as $\lambda_{1}$,
$\lambda_{2}$($\lambda_{1}>\lambda_{2}$). We define local linearity $\theta$
as

$\theta=\frac{\lambda_{1}-\lambda_{2}}{\lambda_{1}}.$ (2)

Lastly, a lower bound $\theta_{min}$ is set to select features with high local
linearity. In order to extract features with high aggregation degree, we
compute the radius(marked as r) of $\mathcal{N}_{i}$, which is the maximum
distance between the feature and the point in $\mathcal{N}_{i}$. If a
feature’s local linearity $\theta>\theta_{min}$ and radius r is less than the
maximum distance $\textbf{d}_{max}$, we consider it as a surface feature point
and keep it. We set $\textbf{m}=10$, $\textbf{d}_{max}=2$, $\theta_{min}=0.9$.

###  III-C Pose Tracking

####  III-C1 Motion Compensation

Since the radar feature points in a scan are received at different times, we
need to align all feature points to the beginning time of the scan. The
process is called Motion Compensation. As an approximation, the radar motion
is modeled with constant angular and linear velocities between two scans. This
allows us to use linear interpolation to calculate the positions of features
at the beginning of the scan.

####  III-C2 Motion Estimation

The pose estimation aligns the current surface feature point cloud
$\mathcal{S}_{k}$ with the global feature map $\mathcal{M}_{k-1}$. We solve
the pose (marked as $T_{k}$) of the kth scan by optimizing the loss function
of the distance between surface feature point cloud $\mathcal{S}_{k}$ and the
global feature map $\mathcal{M}_{k-1}$, and the distance is calculated by
point-to-line metric:

$\displaystyle T_{k}$ $\displaystyle=\arg\min_{T}f(\mathcal{M}_{k-1},\
\mathcal{S}_{k},\ T)$
$\displaystyle=\arg\min_{T}\sum_{\textbf{p}_{i}^{S}\in\mathcal{S}_{k}}||((T\textbf{p}_{i}^{S}-\textbf{p}_{i}^{M})\times\vec{n}_{i}||^{2},$
(3)

where $\textbf{p}_{i}^{S}$ is a surface feature in $\mathcal{S}_{k}$,
$\textbf{p}_{i}^{M}$ is the geometric center of the corresponding points set
in global feature map $\mathcal{M}_{k-1}$, and $\vec{n}_{i}$ is the unit
vector of the norm of the corresponding points set.

###  III-D Probability-based Feature Filter

The feature points selected by geometry-based feature selection still contain
some noise caused by multi-path effect and power saturation. These noise has
little effect when using scan-to-scan or scan-to-frames matching strategy.
However, when _MAROAM_ need to build a local map and use scan-to-map matching
strategy, wrong features will be accumulated in the local map, which will
significantly reduce the accuracy of motion estimation.

Inspired by _pFilter_ [25], we design a probability-based feature filter to
further extract more stable feature points. _pFilter_ propose a new metric
”p-Index” to evaluate feature points, which indicates whether a point is
frequent or transient. With the same thought, we define a simple metric as the
frequency of a feature point being hit.

For a feature point $p_{k_{0}}^{M}$ that appears for the first time in the
$k_{0}$th scan map, we define $\textbf{I}(\textbf{p}_{k_{0}}^{M},k)$ to
indicate whether $\textbf{p}_{k_{0}}^{M}$ is hit in the $k$th ($k>k_{0}$)
scan:

$\textbf{I}(\textbf{p}_{k_{0}}^{M},k)=\left\\{\begin{array}[]{ll}1&\exists\,\textbf{p}_{k}^{S}\in\mathcal{S}_{k},s.t.\
\textbf{p}_{k_{0}}^{M}\in\mathcal{N}_{k}(\textbf{p}_{k}^{S}),\\\
0&\text{otherwise},\end{array}\right.$ (4)

where $\mathcal{S}_{k}$ denotes point cloud of extracted surface feature
points from the $k$th images of radar, and
$\mathcal{N}_{k}(\textbf{p}_{k}^{S})$ denotes the set of corresponding feature
points in $\mathcal{M}_{k}$ for $\textbf{p}_{k}^{S}$ in $S_{k}$. Intuitively,
$\textbf{I}(\textbf{p}_{k_{0}}^{M},k)=1$ indicates that there exits an
extracted feature point in the $k$th scan that can be matched with the point
$\textbf{p}_{k_{0}}^{M}$, and it is called to be hit.

For each feature point $\textbf{p}_{k_{0}}^{M}$ in the feature map
$\mathcal{M}_{k}$, we define $R(\textbf{p}_{k_{0}}^{M},k)$ and
$H(\textbf{p}_{k_{0}}^{M},k)$ to describe the attributes of feature points.
$R(\textbf{p}_{k_{0}}^{M},k)$ represents the matching rounds of the feature
point:

$R(\textbf{p}_{k_{0}}^{M},k)=k-k_{0}+1.$ (5)

And $H(\textbf{p}_{k_{0}}^{M},k)$ represents the total number of hits of the
feature in all historical matches:

$H(\textbf{p}_{k_{0}}^{M},k)=\sum_{k^{\prime}=k_{0}+1}^{k}\textbf{I}(\textbf{p}_{k_{0}}^{M},k^{\prime}).$
(6)

The above formulas can be expressed in recursive form:

$\begin{array}[]{ll}R(\textbf{p}_{k_{0}}^{M},k+1)&=R(\textbf{p}_{k_{0}}^{M},k)+1,\\\
H(\textbf{p}_{k_{0}}^{M},k+1)&=H(\textbf{p}_{k_{0}}^{M},k)+\textbf{I}(\textbf{p}_{k_{0}}^{M},k+1).\end{array}$
(7)

At last, we define the hit probability $P(\textbf{p}_{k_{0}}^{M},k)$ for each
feature $\textbf{p}_{k_{0}}^{M}$ as the hit number divide the round number:

$P(\textbf{p}_{k_{0}}^{M},k)=H(\textbf{p}_{k_{0}}^{M},k)/R(\textbf{p}_{k_{0}}^{M},k).$
(8)

For feature points with matching map feature points in the neighborhood, we
estimate an initial value of $R(\textbf{p}_{k_{0}}^{M},k_{0})$ and
$H(\textbf{p}_{k_{0}}^{M},k_{0})$ for each feature point
$\textbf{p}_{k_{0}}^{M}$, using the average value of $R$ and $H$ for five
feature points that correspond to $\textbf{p}_{k}^{M}$ in the match.

We define three parameters $\theta_{p}$, $R_{min}$, $H_{max}$ for feature
selection:

  * •

If $\textbf{p}_{k_{0}}^{M}$ has just been detected recently, i.e.,
$R(\textbf{p}_{k_{0}}^{M},k)<R_{min}$, then we maintain it in the local map
for a while to see whether it would be stable;

  * •

If a feature $\textbf{p}_{k_{0}}^{M}$ is detected many times, i.e.,
$H(\textbf{p}_{k_{0}}^{M},k)>H_{max}$, we still consider it a stable feature
even if it is not detected for many scans afterward;

  * •

At last, $\theta_{p}$ represents as the threshold to identify stable feature
points. For a feature $\textbf{p}_{k_{0}}^{M}$, If
$H(\textbf{p}_{k_{0}}^{M},k_{0})>\theta_{p}$, it will be maintained in the
local feature map permanently. Otherwise, it will be removed from local
feature map.

The algorithm will select features in two stages: feature matching and map
update. The process of probability-based feature selection is summarized in
Algorithm 1.

Data: $\mathcal{M}_{k},\ \mathcal{S}_{k},\ \textbf{T}_{k}$

Result: $\mathcal{M}_{k+1}$

1 for  _$\textbf{p}_{k}^{S}\emph{{in}}\mathcal{S}_{k}$_ do

2 $\textbf{p}_{k}^{S}\leftarrow\textbf{T}_{k}\textbf{p}_{k}$;

3 $\mathcal{N}_{k}(\textbf{p}_{k})\leftarrow
Find\\_correspondence(\textbf{p}_{k}^{S},\mathcal{M}_{k},5$);

4 for
_$\textbf{p}_{k_{0}}^{M}\textbf{in}\mathcal{N}_{k}(\textbf{p}_{k}^{S})$_ do

5 $H(\textbf{p}_{k_{0}}^{M},k)\leftarrow H(\textbf{p}_{k_{0}}^{M},k)+1$;

6

7 end for

8 $H(\textbf{p}_{k}^{S},k)\leftarrow Average(H(\textbf{p}_{k_{0}}^{M},k)$ for
$\textbf{p}_{k_{0}}^{M}\textbf{in}\mathcal{N}_{k}(\textbf{p}_{k}^{S}))$;

9 $R(\textbf{p}_{k}^{S},k)\leftarrow Average(R(\textbf{p}_{k_{0}}^{M},k)$ for
$\textbf{p}_{k_{0}}^{M}\textbf{in}\mathcal{N}_{k}(\textbf{p}_{k}^{S}))$;

10

11 end for

12$\mathcal{M}_{k+1}\leftarrow\mathcal{M}_{k}\cup\mathcal{S}_{k}$ ;

13 for  _$\textbf{p}_{k_{0}}^{M}\emph{{in}}\mathcal{M}_{k+1}$_ do

14 $R(\textbf{p}_{k_{0}}^{M},k+1)\leftarrow R(\textbf{p}_{k_{0}}^{M},k)+1$;

15 $P(\textbf{p}_{k_{0}}^{M},k+1)\leftarrow
H(\textbf{p}_{k_{0}}^{M},k+1)/R(\textbf{p}_{k_{0}}^{M},k+1)$;

16 if  _$P(\textbf{p}_{k_{0}}^{M},k+1) <\theta_{p}\ \textbf{and}\
R(\textbf{p}_{k_{0}}^{M},k+1)>R_{min}\ \textbf{and}\
H(\textbf{p}_{k_{0}}^{M},k+1)<H_{max}$_ then

17 Delete($\textbf{p}_{k_{0}}^{M}$);

18

19

20 end for

Algorithm 1 Probability-based Feature Filter

##  IV Experimental and Evaluation

We firstly evaluate the performance of _MAROAM_ on the Oxford Radar RobotCar
Dataset[7] which has urban sceneries and variety of weather conditions. Then
we evaluate _MAROAM_ on the MulRan[15] Dataset and the Boreas Dataset[16]
which have different Radar model, non-urban sceneries and challenging weather
conditions. Then, we construct an ablation study to examine the effects of
geometry-probability two-step feature filter. To this end, we also implemented
the scan-to-frames methods to prove that our map-based method has less drift.
All experiments are conducted on a laptop with an AMD Ryzen 9 4800H CPU and
16GB memory.

TABLE I: Evaluation on Oxford Radar RobotCar dataset

Method | 10-12-32 | 16-13-09 | 17-13-26 | 18-14-14 | 18-15-20 | 10-11-46 | 16-11-53 | 18-14-46 | Mean
---|---|---|---|---|---|---|---|---|---
SuMa[26](Lidar) | 1.1/0.3* | 1.2/0.4* | 1.1/0.3* | 0.9/0.1* | 1.0/0.2* | 1.1/0.3* | 0.9/0.3* | 1.0/0.1* | 1.16/0.3*
ORB_SLAM2(Vision) | 6.09/1.6 | 6.23/1.7 | 6.41/1.7 | 7.05/1.8 | 11.5/3.3 | 6.11/1.7 | 6.16/1.7 | 7.17/1.9 | 7.09/3.1
Cen[3] | N/A | N/A | N/A | N/A | N/A | N/A | N/A | N/A | 3.72/0.95
Under the Radar[4] | N/A | N/A | N/A | N/A | N/A | N/A | N/A | N/A | 2.05/0.67**
Mask by Moving[5] | N/A | N/A | N/A | N/A | N/A | N/A | N/A | N/A | 2.78/0.85
HERO[23] | 1.77/0.62 | 1.75/0.59 | 2.04/0.73 | 1.83/0.61 | 2.20/0.77 | 2.14/0.71 | 2.01/0.61 | 1.97/0.65 | 1.96/0.66
Hong odometry[20] | 2.98/0.8 | 3.12/0.9 | 2.92/0.8 | 3.18/0.9 | 2.85/0.9 | 3.26/0.9 | 3.28/0.9 | 3.33/1 | 3.11/0.9
Hong SLAM[20] | 2.17/0.67 | 1.84/0.59 | 2.46/0.81 | 2.21/0.71 | 2.45/0.78 | 2.27/0.9 | 2.24/0.6 | 2.34/0.7 | 2.21/0.7
Hong odometry V2[21] | 2.32/0.7 | 2.62/0.7 | 2.27/0.6 | 2.29/0.7 | 2.25/0.7 | 2.16/0.6 | 2.49/0.7 | 2.12/0.6 | 2.32/0.7
Hong SLAM V2[21] | 1.98/0.6 | 1.48/0.5 | 1.71/0.5 | 2.22/0.7 | 1.77/0.6 | 1.96/0.7 | 1.81/0.6 | 1.68/0.5 | 1.83/0.6
CFEAR[10] | 1.64/0.48 | 1.86/0.52 | 1.66/0.48 | 1.71/0.49 | 1.75/0.51 | 1.65/0.48 | 1.99/0.53 | 1.79/0.50 | 1.76/0.50
_MAROAM_ (w/o loop) | 1.80/0.49 |  1.83/0.57 | 1.92/0.49 | 1.61/0.49 | 1.69/0.51 | 1.73/0.52 | 1.81/0.58 | 1.92/0.53 | 1.79/0.52
_MAROAM_ | 1.63/0.46 |  1.83/0.56 | 1.49/0.47 | 1.54/0.47 | 1.61/0.50 |  1.55/0.53 |  1.78/0.54 | 1.55/0.50 | 1.62/0.50

  * •

TBALE I shows the evaluation results of different algorithms on Oxford Radar
RobotCar Dataset[7]. Bold font indicate the best results. Mark N/A means the
algorithm cannot finish the full sequence. Results marked * cannot be compared
directly as these can’t run the whole trajectory and fail at some time in the
middle. Results marked ** cannot be compared directly as these are trained and
evaluated on the same spatial location.

###  IV-A Performance on Oxford Radar RobotCar Dataset

The Oxford Radar RobotCar Dataset[7] is an open large-scale radar dataset used
for evaluating Radar-based SLAM. It includes 32 sequences of radar data
collected while traversing a same route in Oxford with ground truth. The radar
data is captured by a Navtech CTS-350X, a Frequency Modulated Continuous Wave
(FMCW) scanning radar, which is configured to return 3768 power readings at a
resolution of 4.32cm across 400 azimuths and operates at the frequency of 4Hz.
For comparison, we selected the same 8 sequences as CFEAR Radarodometry[10]
and RadarSLAM[20]. For all sequences, we set $\theta_{m}=0.25$, $R_{min}=10$
and $H_{max}=10$. The odometry accuracy was measured as proposed in the KITTI
odometry benchmark[27] to compute the average translation error (%) and
rotation error (deg/m) over all sub sequences between 100, 200, . . . , 800 m.

We compared _MAROAM_ with 7 state-of-the-art radar based odometry and SLAM
algorithms. The evaluation result is shown in TABLE I. Experiment shows that
_MAROAM_ achieves the best accuracy with a translation of 1.62% and rotation
error of 0.50 deg/100m. The accuracy of _MAROAM_ (1.62%) is better than
CFEAR[10] (1.76%), RadarSLAM[21, 20] (1.83%, 2.21%), HERO[23] (1.96%), Mask by
Moving[5] (without CSV, 2.78%), Under the Radar[4] (2.05%). Without loop
closure, the mean accuracy of our odometry (1.79%) is close to that of the
SOTA algorithm CFEAR (1.76%). When using loop closure, _MAROAM_ imporve by
7.95% relative to CFEAR.

###  IV-B Performance on MulRan Dataset and Boreas Dataset

TABLE II: ATE trans. (RMSE) [m] on MulRan Dataset Method | Evaluation | DCC02 | Riverside02 | Mean
---|---|---|---|---
PhaRaO[18] |  | 24.962 | 31.83 | 28.40
SuMa[26] | [21] | 17.834 | N/A | N/A
Hong SLAM [20] |  | 24.962 | 95.247 | 60.10
Hong SLAM v2 [21] |  | 9.878 | 7.049 | 8.46
_MAROAM_ (w/o loop) |  | 7.52 | 16.33 | 11.93
_MAROAM_ |  | 5.81 | 4.85 | 5.33
TABLE III: EVALUATION ON BOREAS DATASET Algorithm | Sequence | Weather | Error | Mean Error
---|---|---|---|---
HERO[23] | 01-26-10-59 | snow | 2.00/0.56 | 2.02/0.56
01-26-11-22 | snow | 1.98/0.53
02-09-12-55 | sun | 2.07/0.59
_MAROAM_ (w/o loop) | 11-26-13-58 | snow | 1.78/0.44 | 1.84/0.45
12-01-13-26 | snow | 2.02/0.45
12-18-13-44 | sun | 1.72/0.45
Figure 4: Estimated Odometry Trajectories, SLAM Trajectories and Ground Truth
of 2 Sequences from MulRan Dataset[15].

(a) 11-26-13-58

(b) 12-01-13-26

(c) 12-18-13-44

Figure 5: Estimated odometry (without loop closure) trajectories and
groundtruth of 3 sequences from Boreas dataset[15].

To test the performance of _MAROAM_ in non-urban scenery and challenging
weathers, we evaluate it on MulRan Dataset[15] and Boreas Dataset[16]. The
radar images in both datasets are collected from Navtech CIR204-H FMCW
scanning radar which is different from Oxford dataset used. This Navtech
CIR204-H radar has a range resolution of 0.0596 m, and a total range of 200 m,
which is different from Navtech CTS-350X. For both datasets, we set
$\theta_{m}=0.2$, $R_{min}=10$ and $H_{max}=10$.

####  IV-B1 MulRan

the MulRan dataset[15] dataset includes four different scenarios of urbans and
non-suburbs. In order to compare with PhaRao[18] and Hong SLAM[20, 21], we
selected DCC02 and Riverside02 sequences for evaluation. The third sequence in
PhaRao is not open sourced in MulRan. We use the open source tool evo[28] to
evaluate the absolute trajectory error (ATE) of the trajectory.

The trajectory and quantitative results on the sequences are shown in Fig. 4
and TABLE II. Experiments show that _MAROAM_ achieves the best accuracy with
ATE of 5.33 m. For comparing the accuracy of the odometry (w/o loop), _MAROAM_
(w/o loop) (11.93m) is much better than that of PhaPaO (28.4 m). Compared with
the complete SLAM pipeline, the drift error of _MAROAM_ (5.33 m) also exceeds
that of Hong SLAM (60.10 m) and Hong SLAM v2 (8.46 m). _MAROAM_ reduces drift
error by 37.0% relative to the SOTA algorithm Hong SLAM v2. Among them, in
Riverside02 sequence, which is open, empty and lacks structural textures, our
odometry and SLAM algorithms can operate stably where SuMa[26] can’t finish
the full sequence in Hong[21]’s evaluation.

####  IV-B2 Boreas

The Boreas Dataset[16] is the latest open source dataset with diverse weather
conditions for SLAM and object detection research. In order to compare with
HREO[23], we selected three sequences for evaluation. Like HERO’s selction,
two sequences were taken during a snow day, and the other was taken on a sunny
day. The odometry accuracy was measured as proposed in the KITTI odometry
benchmark[27]. To be fair, we only tested _MAROAM_ with loop closure disabled.

The trajectory and quantitative results on the sequences are shown in Fig. 5
and TABLE III. Experiment shows that _MAROAM_ achieves the accuracy with a
translation of 1.84% and rotation error of 0.45 deg/100m, exceeds HREO’s 2.02%
and 0.56 deg/100m. _MAROAM_ reduces translation error by 8.9% and rotation
error by 19.6% relative to HERO.

###  IV-C Ablation Study

Figure 6: Trajectories of ablation study using sequence ”10-12-32” of the Oxford dataset.. TABLE IV: Ablation Experiment | Trans.(%) | Rot.(rad/100m)
---|---|---
scan-to-frames(5) | 2.9 | 0.84
scan-to-frames(7) | 2.57 | 0.74
scan-to-frames(9) | 2.52 | 0.73
_MAROAM_ (w/o probability filter) | 3.16 | 0.97
_MAROAM_ (w/o geometry filter) | 8.0 | 2.86
_MAROAM_ (w/o loop) | 1.80 | 0.49
_MAROAM_ (w loop) | 1.63 | 0.46

We show how the design of each module in the system affects the performance of
the proposed framework using sequence ”10-12-32” of the Oxford dataset .The
experimental results are shown in TABLE LABEL:table2 and Fig. 6.

####  IV-C1 Effect of map-based method

We first compare our map-based radar odometry with the scan-to-frames method
which using the sum of the feature points from previous frames as the local
map. The table and figure item scan-to-frames(n) means the algorithm using $n$
frames as the local map which composed of the surface feature points.
Experiments show that the accuracy of our map-based method (1.80%, 0.49) is
28.6% higher than that of scan-to-frames (2.52%, 0.74).

####  IV-C2 Effect of geometry filter

We disable the geometry filter and use the raw feature points to calculate
relative pose and build local map and the odometry error increased from 1.80%,
0.49 rad/100m to 8.00%, 2.86 rad/100m.

####  IV-C3 Effect of probability filter

We disable the probability filter and use the sum of the feature points from
all previous frames as the local map. In our experiments, the odometry error
increased from 1.80%, 0.49 rad/100m to 3.16%, 0.97 rad/100m

####  IV-C4 Effect of loop closure

We have done experiments to observe the effect of loop closure on odometry
accuracy. Experiment shows that the loop closure has obvious effect on
reducing the trajectory drift error. With loop closure, the error dropped from
1.80%, 0.49 rad/100m to 1.63%, 0.46 rad/100m.

##  V Conclusion

In this letter, we propose _MAROAM_ , a framework for map-based radar odometry
and mapping for performing state estimation and mapping in complex scenery,
weather, and road conditions. The proposed framework uses geometry-probability
two-step feature selection in the registration and map update stages. The
algorithm first extract surface feature points by calculating the local
curvature and local aggregation. A point-to-line ICP method then uses surface
feature points to solve the relative pose transformation. Meanwhile, a
probability-based filter is used to dynamically filter out the feature points
that are not hit frequently in the ICP matching process. We evaluate our
_MAROAM_ on three datasets. Experiments show that our algorithm outperforms
SOTA algorithm on all three datasets. As devoted contributors to the open-
source community, we will open source the algorithm after the paper is
accepted.

## References

  * [1] S. Cen, “Ego-motion estimation and localization with millimeter-wave scanning radar,” Master’s thesis, University of Oxford, 2020.
  * [2] S. H. Cen and P. Newman, “Precise ego-motion estimation with millimeter-wave radar under diverse and challenging conditions,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , 2018, pp. 6045–6052.
  * [3] ——, “Radar-only ego-motion estimation in difficult settings via graph matching,” in _2019 International Conference on Robotics and Automation (ICRA)_. IEEE, 2019, pp. 298–304.
  * [4] D. Barnes and I. Posner, “Under the radar: Learning to predict robust keypoints for odometry estimation and metric localisation in radar,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_ , 2020, pp. 9484–9490.
  * [5] D. Barnes, R. Weston, and I. Posner, “Masking by moving: Learning distraction-free radar odometry from pose information,” 2019. [Online]. Available: https://arxiv.org/abs/1909.03752
  * [6] C. Cadena, L. Carlone, H. Carrillo, Y. Latif, D. Scaramuzza, J. Neira, I. Reid, and J. J. Leonard, “Past, present, and future of simultaneous localization and mapping: Toward the robust-perception age,” _IEEE Transactions on Robotics_ , vol. 32, no. 6, pp. 1309–1332, 2016.
  * [7] D. Barnes, M. Gadd, P. Murcutt, P. Newman, and I. Posner, “The oxford radar robotcar dataset: A radar extension to the oxford robotcar dataset,” in _Proceedings of the IEEE International Conference on Robotics and Automation (ICRA)_ , Paris, 2020. [Online]. Available: https://arxiv.org/abs/1909.01300
  * [8] V. Borkar, A. Ghosh, R. Singh, and N. Chourasia, “Radar cross-section measurement techniques,” _Defence Science Journal_ , vol. 60, no. 2, p. 204, 2010.
  * [9] P. J. Besl and N. D. McKay, “A method for registration of 3-d shapes,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 14, pp. 239–256, 1992.
  * [10] D. Adolfsson, M. Magnusson, A. Alhashimi, A. J. Lilienthal, and H. Andreasson, “Cfear radarodometry – conservative filtering for efficient and accurate radar odometry,” 2021.
  * [11] P.-C. Kung, C.-C. Wang, and W.-C. Lin, “A normal distribution transform-based radar odometry designed for scanning and automotive radars,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , 2021, pp. 14 417–14 423.
  * [12] X. Huang, G. Mei, J. Zhang, and R. Abbas, “A comprehensive survey on point cloud registration,” 2021. [Online]. Available: https://arxiv.org/abs/2103.02690
  * [13] J. Zhang and S. Singh, “Loam: Lidar odometry and mapping in real-time,” in _Robotics: Science and Systems_ , 2014.
  * [14] G. Kim and A. Kim, “Scan context: Egocentric spatial descriptor for place recognition within 3d point cloud map,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2018, pp. 4802–4809.
  * [15] G. Kim, Y. S. Park, Y. Cho, J. Jeong, and A. Kim, “Mulran: Multimodal range dataset for urban place recognition,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_ , 2020, pp. 6246–6253.
  * [16] K. Burnett, D. J. Yoon, Y. Wu, A. Z. Li, H. Zhang, S. Lu, J. Qian, W.-K. Tseng, A. Lambert, K. Y. Leung, A. P. Schoellig, and T. D. Barfoot, “Boreas: A multi-season autonomous driving dataset,” _arXiv preprint arXiv:2203.10168_ , 2022.
  * [17] P. Checchin, F. Gérossier, C. Blanc, R. Chapuis, and L. Trassoudaine, “Radar scan matching slam using the fourier-mellin transform,” in _Field and Service Robotics_ , A. Howard, K. Iagnemma, and A. Kelly, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2010, pp. 151–161.
  * [18] Y. S. Park, Y.-S. Shin, and A. Kim, “Pharao: Direct radar odometry using phase correlation,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_ , 2020, pp. 2617–2623.
  * [19] P.-C. Kung, C.-C. Wang, and W.-C. Lin, “A normal distribution transform-based radar odometry designed for scanning and automotive radars,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , 2021, pp. 14 417–14 423.
  * [20] Z. Hong, Y. R. Pétillot, A. M. Wallace, and S. Wang, “Radar slam: A robust slam system for all weather conditions,” _ArXiv_ , vol. abs/2104.05347, 2021\.
  * [21] Z. Hong, Y. R. Petillot, A. M. Wallace, and S. Wang, “Radar SLAM: A robust SLAM system for all weather conditions,” _CoRR_ , vol. abs/2104.05347, 2021\. [Online]. Available: https://arxiv.org/abs/2104.05347
  * [22] L. He, X. Wang, and H. Zhang, “M2dp: A novel 3d point cloud descriptor and its application in loop closure detection,” in _2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , Oct 2016, pp. 231–237.
  * [23] K. Burnett, D. J. Yoon, A. P. Schoellig, and T. D. Barfoot, “Radar odometry combining probabilistic estimation and unsupervised feature learning,” 2021. [Online]. Available: https://arxiv.org/abs/2105.14152
  * [24] H. Wang, C. Wang, C. Chen, and L. Xie, “F-loam : Fast lidar odometry and mapping,” in _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2020.
  * [25] Y. Duan, J. Peng, Y. Zhang, J. Ji, and Y. Zhang, “Pfilter: Building persistent maps through feature filtering for fast and accurate lidar-based slam,” 2022\. [Online]. Available: https://arxiv.org/abs/2208.14848
  * [26] X. Chen, A. Milioto, E. Palazzolo, P. Giguère, J. Behley, and C. Stachniss, “Suma++: Efficient lidar-based semantic slam,” in _2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2019, pp. 4530–4537.
  * [27] A. Geiger, P. Lenz, and R. Urtasun, “Are we ready for autonomous driving? the kitti vision benchmark suite,” in _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2012.
  * [28] M. Grupp, “evo: Python package for the evaluation of odometry and slam.” https://github.com/MichaelGrupp/evo, 2017.


# Quadratic Optimization-Based Nonlinear Control for Protein Conformation
Prediction

Alireza Mohammadi and Mark W. Spong A. Mohammadi is with the Department of
Electrical & Computer Engineering, University of Michigan-Dearborn, Dearborn,
MI 48128 USA. Mark W. Spong is with the University of Texas at Dallas,
Richardson, TX 75080 USA. Mark W. Spong is the Excellence in Education Chair
at the University of Texas at Dallas. Emails<EMAIL_ADDRESS><EMAIL_ADDRESS>∗Corresponding Author: A. Mohammadi

###### Abstract

Predicting the final folded structure of protein molecules and simulating
their folding pathways is of crucial importance for designing viral drugs and
studying diseases such as Alzheimer’s at the molecular level. To this end,
this paper investigates the problem of protein conformation prediction under
the constraint of avoiding high-entropy-loss routes during folding. Using the
well-established kinetostatic compliance (KCM)-based nonlinear dynamics of a
protein molecule, this paper formulates the protein conformation prediction as
a pointwise optimal control synthesis problem cast as a quadratic program
(QP). It is shown that the KCM torques in the protein folding literature can
be utilized for defining a reference vector field for the QP-based control
generation problem. The resulting kinetostatic control torque inputs will be
close to the KCM-based reference vector field and guaranteed to be constrained
by a predetermined bound; hence, high-entropy-loss routes during folding are
avoided while the energy of the molecule is decreased.

##  I Introduction

Computer-aided prediction of the folded structure of a protein molecule lies
at the heart of protein engineering, drug discovery, and investigating
diseases such as Alzheimer’s at the molecular and cellular levels [1]. The 3D
structure of a protein molecule, known as the protein conformation, is mainly
determined by its linear amino acid (AA) sequence [2]. The protein folding
problem is concerned with determining the final folded structure of a protein
molecule given its linear AA sequence. Studying the folding pathways is also
important for designing viral drugs that cause misfolding in virus proteins
[3].

Knowledge-based and physics-based methods are the two prominent computational
approaches for solving the protein folding problem. Knowledge-based methods,
which also include the approach of Google’s DeepMind AlphaFold [4], rely on
using previously determined types of folds to solve the protein folding
problem [5]. In physics-based methods, on the other hand, protein folding
simulations are carried out using the first principles [6]. Physics-based
methods enjoy several advantages over their knowledge-based counterparts such
as the ability to model the protein-nucleic acid interactions and to _simulate
protein folding pathways_. However, physics-based methods relying on standard
molecular dynamics (MD) suffer from numerical instabilities [7]. To address
the inherent challenges of MD-based protein folding solutions, Kazerounian and
collaborators [8, 9, 10, 11] introduced the kinetostatic compliance method
(KCM) to overcome the computational problems arising from large degrees-of-
freedom (DOFs) and to provide a faster convergence to the molecular minimum
energy state. This framework assumes that folding takes place as a quasi-
equilibrium process where the protein chain complies under the kinetostatic
effect of the inter-atomic and intramolecular force fields. The KCM framework
has been successful in investigating phenomenon such as the effect of hydrogen
bond formation on the protein mechanical mobility [12] and conformational
trajectory planning of proteins in 3D workspace [13].

Along various modeling approaches, numerous optimization techniques have also
been adopted for obtaining realistic folding pathways such as homotopy-based
[14] and optimal control-based [15, 16, 17] algorithms. Using a bead-spring
linear time-invariant (LTI) model of the protein molecule, Arkun and Erman in
[15] propose using a linear quadratic regulator (LQR) control input for
guiding the folding simulations in order to simultaneously evade high-energy
regions of the folding landscape and to _avoid_ high-entropy-loss routes. An
interesting aspect of Arkun and Erman’s optimal control framework, which is
inspired from the behavior of the two-state protein molecules [18], is to
provide _the optimal solution to protein folding as a trade-off between the
molecule energy minimization and the avoidance of high-entropy-loss routes_.
In particular, to avoid extremely fast decay of excess entropy, Arkun and
Erman penalize the control input by using an infinite horizon LQR cost
function. Arkun and Gür in [16] extended the aforementioned work by providing
the sampled trajectories of the coarse-grained protein model as initial
conditions for an all-atom MD-based simulation environment, namely, Nanoscale
Molecular Dynamics (NAMD) software [19].

Our point of departure in this paper is a synergy of nonlinear KCM-based
folding dynamics due to Kazerounian and collaborators [9, 10] and LQR control
of folding due to Arkun and collaborators for a linear time-invariant system
[15, 16]. In particular, we propose using a pointwise optimal control law for
guiding the protein folding simulations under the nonlinear KCM assumptions.
Unlike the bead-spring LTI model in [15, 16] that uses the Cartesian
coordinates of the atoms, the KCM model in [9, 10] directly accounts for the
effect of the nonlinear interatomic interactions on the protein molecule
dihedral angles. Our feedback control solution is based on the quadratic
optimization of control torques, by using the optimal decision strategy (ODS)
framework due to Barnard [20]. Our QP-based control inputs minimize the
deviation between the open-loop dynamics vector field and a reference model
vector field under the entropy-loss constraints during the folding procedure.
Optimal decision strategies have been employed in applications ranging from
the control of manipulators with bounded inputs [21] and magnetic microrobots
under frequency constraints [22] to controlling transients in power systems
[23]. ODS-based strategies belong to the larger family of optimization-based
nonlinear controllers [24, 25, 26, 27], whose applications in robotics are
growing, thanks in part to advancements in mobile computation power.
Contributions of the paper. This paper contributes to solving the protein
conformation prediction and obtaining more realistic folding pathways using
the KCM-based framework in two important ways. First, we demonstrate how the
KCM-based protein folding difference equation can be cast as a control
synthesis problem and that the kinetostatic torque in [9, 10] can be
interpreted as a control input that guides the folding simulation (Proposition
III.2). Second, using the control torque derived from the work in [9, 10], we
define a KCM-based reference vector field for the folding dynamics. This KCM-
based vector field will then be used in a quadratic program (QP) for
synthesizing bounded control inputs with predetermined bounds. Therefore, our
synthesized control inputs will not only be close to the KCM-based reference
vector field in [9, 10] but also bounded by user-specified bounds, thus
respecting the requirement of avoiding high-entropy-loss routes during the
folding process. Finally, Proposition IV.1 provides sufficient conditions for
uniqueness and Lipschitz continuity of the control inputs generated by our
proposed QP.

The rest of this paper is organized as follows. First, we briefly review the
kinematics and dynamics of protein folding in Section II. Next, we demonstrate
how the KCM-based approach to protein folding can be formulated as a control
synthesis problem in Section III. Thereafter, we present our ODS-based control
scheme for protein folding in Section IV. After presenting the simulation
results in Section V, we conclude the paper with final remarks and future
research directions in Section VI.

Notation. We let $\mathbb{R}_{+}$ denote the set of all non-negative real
numbers. Given a real number $x$, we let $\lfloor x\rfloor$ denote the
greatest integer less than or equal to $x$. Given
$\mathbf{x}\in\mathbb{R}^{M}$, we let
$|\mathbf{x}|:=\sqrt{\mathbf{x}^{\top}\mathbf{x}}$ and
$|\mathbf{x}|_{\infty}:=\max\limits_{i}|x_{i}|$ denote the Euclidean and the
maximum norms of $\mathbf{x}$, respectively. Given the integer $M$, we let
$\mathbf{e}_{M}$ and $\mathbf{I}_{M}$ denote the vector
$[1,\,\cdots,\,1]^{\top}$ and the identity matrix of size $M$, respectively.

##  II Preliminaries

In this section, we briefly review the KCM framework for protein folding _in
vacuo_. For brevity, our presentation will be focused on the protein main
chain. This simplification is supported by the observation that the motion of
the backbone $\text{C}_{\alpha}$ atoms explains most of the essential folding
dynamics (see, e.g., [28]). However, certain phenomenon such as hydrophobic
core collapse can only be explained by considering the side chains [29].
Details on the inclusion of side chains can be found in [10, 9, 8].

###  II-A Kinematic linkage representation of protein molecules

Proteins are long molecular chains consisting of peptide planes that are
rigidly joined together via peptide bonds (see Figure 1). Except the first and
the last peptide planes, each pair of these planes are hinged to each other
via a central carbon atom, known as the alpha-Carbon atom and denoted by
$\text{C}_{\alpha}$. The coplanarity assumption of the six atoms
$\text{C}_{\alpha}-\text{CO}-\text{NH}-\text{C}_{\alpha}$ in each amino acid
is due to the experimental observations of the structure of protein molecules
using high resolution X-ray crystallographic techniques [2]. In Figure 1, the
red line segments represent covalent chemical bonds between the atoms. A
hydrogen atom is also covalently bonded to each $\text{C}_{\alpha}$ atom,
which is not shown in the figure.

Figure 1: The protein molecule kinematic structure under the
$\text{C}_{\alpha}-\text{CO}-\text{NH}-\text{C}_{\alpha}$ coplanarity
assumption. For a large-scaled version of this figure, see Supplementary
Material.

As it can be seen from Figure 1, each $\text{C}_{\alpha}$ is connected to four
other chemical components, namely, the three atoms N, H, and C as well as a
variable side chain denoted by SR. The first $\text{C}_{\alpha}$ of the
linkage is connected to an amino group, known as the N-terminus, and one other
peptide plane. The last $\text{C}_{\alpha}$ of the linkage hinges to a
carboxyl group, known as the C-terminus, and one other peptide plane.

The backbone conformation consisting of
$-\text{N}-\text{C}_{\alpha}-\text{C}-$ atoms is completely determined using a
set of bond lengths and two sets of dihedral angles, i.e., the rotation angles
around $\text{N}-\text{C}_{\alpha}$ and $\text{C}_{\alpha}-\text{C}$.
Therefore, the vector
$\boldsymbol{\theta}=\big{[}\theta_{1},\cdots,\theta_{2N}]^{\top}\in\mathbb{R}^{2N}$
represents the backbone chain configuration with $N-1$ peptide planes.
Corresponding to each DOF of the protein molecule, we consider a unit vector
along its corresponding rotation axis and denote it by $\mathbf{u}_{j}$,
$1\leq j\leq 2N$. Hence, the vectors $\mathbf{u}_{2i}$ and $\mathbf{u}_{2i+1}$
are the unit vectors along the $\text{C}_{\alpha}-\text{C}$ and
$\text{N}-\text{C}_{\alpha}$ bonds of the $i$-th amino acid, respectively.
Moreover, $\mathbf{u}_{1}$ and $\mathbf{u}_{2N}$ correspond to the unit
vectors of the N\- and C-termini, respectively.

In addition to the unit vectors $\mathbf{u}_{j}$, Kazerounian and
collaborators use the so-called body vectors to complete the spatial
orientation description of the rigid peptide planes. The body vectors are
denoted by $\mathbf{b}_{j}$, $1\leq j\leq 2N$ and they completely determine
the relative position of any two atoms in the peptide planes. In particular,
the relative positions of any two atoms can be described by the linear
combination $k_{1m}\mathbf{b}_{2i}+k_{2m}\mathbf{b}_{2i+1}$, where the
constants $k_{1m}$ and $k_{2m}$, $1\leq m\leq 4$, are the same across all
peptide planes (see [30, 11] for further details). Using the vectors
$\mathbf{u}_{j}$ and $\mathbf{b}_{j}$, one can completely describe the
conformation of a protein molecule using the vector of dihedral angles
$\boldsymbol{\theta}$. In particular, after designating the zero position
conformation with $\boldsymbol{\theta}=\mathbf{0}$ corresponding to the
biological reference position of the chain (see [8]), the transformations

$\mathbf{u}_{j}(\boldsymbol{\theta})=\Xi(\boldsymbol{\theta},\mathbf{u}_{j}^{0})\mathbf{u}_{j}^{0},\,\mathbf{b}_{j}(\boldsymbol{\theta})=\Xi(\boldsymbol{\theta},\mathbf{u}_{j}^{0})\mathbf{b}_{j}^{0},$
(1)

where
$\Xi(\boldsymbol{\theta},\mathbf{u}_{j}^{0})=\prod_{r=1}^{j}R(\theta_{j},\mathbf{u}_{j}^{0})$,
can be used to describe the kinematic configuration of the protein molecule
given the vector $\boldsymbol{\theta}$. In (1), the rotation matrix
$R(\theta_{j},\mathbf{u}_{j}^{0})\in\text{SO}(3)$ represents the rotation
about the unit vector $\mathbf{u}_{j}^{0}$ with angle $\theta_{j}$. Having
obtained the body vectors $\mathbf{b}_{j}(\boldsymbol{\theta})$ from (1) and
under the assumption that the N-terminus atom is fixed at the origin, the
coordinates of the backbone chain atoms in the $k$th-peptide plane can be
found from
$\mathbf{r}_{i}(\boldsymbol{\theta})=\sum_{j=1}^{i}\mathbf{b}_{j}(\boldsymbol{\theta}),\;\;1\leq
i\leq 2N-1$ where the indices $i=2k-1$ and $i=2k$ correspond to the N and
$\text{C}_{\alpha}$ atoms, respectively.

###  II-B KCM-based folding

The KCM framework in [9, 10, 11] relies on the well-established fact that the
inertial forces during the folding of a protein chain are negligible with
respect to their electrostatic and van der Waals interatomic counterparts
(see, e.g., [31, 9, 15, 16]). Instead, the dihedral angles change under the
kinetostatic effect of the interatomic force fields with an amount
proportional to the effective torques on these joints.

To present the KCM framework, let us consider a protein chain with $N_{a}$
atoms and $N-1$ peptide planes with the dihedral angle conformation vector
$\boldsymbol{\theta}\in\mathbb{R}^{2N}$ and use the following notation. Given
any two atoms $a_{i}$, $a_{j}$ in the protein chain, we denote their Cartesian
position vectors by $r_{i}(\boldsymbol{\theta})$,
$r_{j}(\boldsymbol{\theta})$, their distance by
$d_{ij}(\boldsymbol{\theta}):=|r_{i}(\boldsymbol{\theta})-r_{j}(\boldsymbol{\theta})|$,
their electrostatic charges by $q_{i}$, $q_{j}$, their van der Waals radii by
$R_{i}$, $R_{j}$, their van der Waals distance by $D_{ij}=R_{i}+R_{j}$, their
dielectric constant by $\varepsilon_{ij}$, and their depth of potential well
as $\epsilon_{ij}$. Finally, we let $w_{ij}^{\text{elec}}$ and
$w_{ij}^{\text{vdw}}$ represent the weight factors for the electrostatic and
van der Waals forces between the two atoms $a_{i}$ and $a_{j}$, respectively.
All of these parameters are reported in [11] and the references therein.

The successive iteration of the KCM-based folding process is performed until
all of the kinetostatic torques converge to a minimum corresponding to a local
minimum of the aggregated free energy function of the protein chain. In
particular, the aggregated free energy is given by

$\mathcal{G}(\boldsymbol{\theta}):=\mathcal{G}^{\text{elec}}(\boldsymbol{\theta})+\mathcal{G}^{\text{vdw}}(\boldsymbol{\theta}),$
(2)

where

$\mathcal{G}^{\text{elec}}(\boldsymbol{\theta})=\sum_{i=1}^{N_{\text{a}}}\sum_{j\neq
i}\frac{w_{ij}^{\text{elec}}}{4\pi\varepsilon_{ij}}\frac{q_{i}q_{j}}{d_{ij}(\boldsymbol{\theta})},$
(3)

is the potential energy of the molecule due to the electrostatic interactions
between the atoms, and

$\mathcal{G}^{\text{vdw}}(\boldsymbol{\theta})=\sum_{i=1}^{N_{\text{a}}}\sum_{j\neq
i}w_{ij}^{\text{vdw}}\epsilon_{ij}\bigg{[}\big{(}\frac{D_{ij}}{d_{ij}(\boldsymbol{\theta})}\big{)}^{12}-2\big{(}\frac{D_{ij}}{d_{ij}(\boldsymbol{\theta})}\big{)}^{6}\bigg{]}$
(4)

is the potential energy due to the van der Waals interactions. The resultant
Coulombic and Van der Waals forces on each atom $a_{i}$, $1\leq i\leq
N_{\text{a}}$, are given by
$F_{i}^{\text{elec}}(\boldsymbol{\theta})=-\nabla_{\mathbf{r}_{i}}\mathcal{G}^{\text{elec}}$
and
$F_{i}^{\text{vdw}}(\boldsymbol{\theta})=-\nabla_{\mathbf{r}_{i}}\mathcal{G}^{\text{vdw}}$,
respectively.

Computing the resultant forces and torques on each of the $N-1$ peptide
planes, which are considered as the rigid links of the peptide linkage, and
appending them in the generalized force vector
$\mathcal{F}(\boldsymbol{\theta})\in\mathbb{R}^{6N}$, one can use a proper
mapping to relate $\mathcal{F}(\boldsymbol{\theta})$ to the equivalent torque
vector acting on the dihedral angles. In particular, starting from an initial
conformation vector $\boldsymbol{\theta}_{0}$, the KCM-based iteration can be
written using the difference equation

$\boldsymbol{\theta}_{k+1}=\boldsymbol{\theta}_{k}+\frac{h}{|\boldsymbol{\tau}(\boldsymbol{\theta}_{k})|_{\infty}}\boldsymbol{\tau}(\boldsymbol{\theta}_{k}),\;$
(5)

where $k$ is a non-negative integer and $h$ is a positive real constant that
tunes the maximum dihedral angle rotation magnitude in each step. The vector
$\boldsymbol{\tau}(\boldsymbol{\theta}_{k})\in\mathbb{R}^{2N}$, which
represents the overall joint torques due to the interatomic forces in the
protein molecule, is given by

$\boldsymbol{\tau}(\boldsymbol{\theta}_{k})=\mathcal{J}^{\top}(\boldsymbol{\theta}_{k})\mathcal{F}(\boldsymbol{\theta}_{k}).$
(6)

In (6), the matrix
$\mathcal{J}(\boldsymbol{\theta}_{k})\in\mathbb{R}^{6N\times 2N}$ represents
the Jacobian of the protein chain at conformation
$\boldsymbol{\theta}=\boldsymbol{\theta}_{k}$ (see [10, 9] for the calculation
details). The vector $\mathcal{F}(\boldsymbol{\theta}_{k})$ is due to the
torques and forces acting on the peptide planes at
$\boldsymbol{\theta}=\boldsymbol{\theta}_{k}$.

##  III KCM-based Folding Iteration as a Nonlinear Control-Affine System

In this section we formulate the KCM-based folding in (5), (6) as a control
synthesis problem. Such control theoretic point of view can be traced back to
the field of targeted molecular dynamics (TMD) [32] (see also [17, 16]) where
the distance of the protein molecule to a target conformation is constrained
by means of a proportional control input for accelerating the transition to
the final structure.

Consider the following nonlinear control-affine system

$\dot{\boldsymbol{\theta}}=\mathcal{J}^{\top}(\boldsymbol{\theta})\mathcal{F}(\boldsymbol{\theta})+\mathbf{u}_{c},$
(7)

where the Jacobian matrix $\mathcal{J}(\boldsymbol{\theta})$ and the
generalized force vector $\mathcal{F}(\boldsymbol{\theta})$ are the same as in
(6). The nonlinear control-affine system in (7) governs the evolution of the
dihedral angles under the effect of interatomic forces and an additional
control input for constraining the folding process simulation.

###### Remark III.1

In the optimal control synthesis considered by Arkun and collaborators [15,
17, 16], the atomic interactions are modeled by linear spring-like forces
coupled via a proper connectivity matrix. In this work, we are considering a
nonlinear dynamical system (as opposed to the linear dynamics in [15, 17, 16])
in the dihedral angle space of the molecule (as opposed to [15, 17, 16], where
the state vector of the system is the Cartesian position vector of the
$\text{C}_{\alpha}$ atoms). The Supplementary Material provides an in-depth
comparison with [15, 17, 16] (Section B).

The following proposition relates the time trajectories of the nonlinear
control-affine system in (7) and the conformation samples obtained from the
iterations in (5) and (6).

###### Proposition III.2

Consider the KCM-based iteration in (5), (6), the nonlinear control-affine
system in (7), and the state feedback control law

$\mathbf{u}_{c}:=\mathbf{u}_{\text{KCM}}(\boldsymbol{\theta})=\frac{1-|\mathcal{J}^{\top}(\boldsymbol{\theta})\mathcal{F}(\boldsymbol{\theta})|_{\infty}}{|\mathcal{J}^{\top}(\boldsymbol{\theta})\mathcal{F}(\boldsymbol{\theta})|_{\infty}}\mathcal{J}^{\top}(\boldsymbol{\theta})\mathcal{F}(\boldsymbol{\theta}).$
(8)

Given an initial conformation $\boldsymbol{\theta}=\boldsymbol{\theta}_{0}$
and a final time $t^{\ast}>0$, assume that
$\sup\big{|}\mathcal{J}^{\top}(\boldsymbol{\theta})\mathcal{F}(\boldsymbol{\theta})\big{|}\leq\lambda$
for some $\lambda>0$ in a neighborhood $\mathcal{N}_{\boldsymbol{\theta}_{0}}$
of $\boldsymbol{\theta}_{0}$. Moreover, suppose that the Lipschitz condition
$\big{|}\mathcal{J}^{\top}(\boldsymbol{\theta})\mathcal{F}(\boldsymbol{\theta})-\mathcal{J}^{\top}(\boldsymbol{\theta}^{\prime})\mathcal{F}(\boldsymbol{\theta}^{\prime})\big{|}\leq\lambda^{\prime}|\boldsymbol{\theta}^{\prime}-\boldsymbol{\theta}|$
holds for some $\lambda^{\prime}>0$, and all $\boldsymbol{\theta}$,
$\boldsymbol{\theta}^{\prime}$ in the neighborhood
$\mathcal{N}_{\boldsymbol{\theta}_{0}}$. Then, the trajectory of (7) starting
from $\boldsymbol{\theta}_{0}$ under the control law in (8), as long as it
remains in $\mathcal{N}_{\boldsymbol{\theta}_{0}}$ during the time interval
$[0,\,t^{\ast}]$, are related to the solutions of the KCM-based difference
equation (5), (6) through

$\displaystyle\big{|}\boldsymbol{\theta}(hk)-\boldsymbol{\theta}_{k}\big{|}\leq\frac{1+\lambda}{2}(\exp(t^{\ast}\lambda^{\prime})-1)h,\;$
$\displaystyle\text{ for all
}k=0,\,1,\,\cdots,\,\lfloor\frac{t^{\ast}}{h}\rfloor$ (9)

###### Proof:

See Section D in Supplementary Material. ∎

###### Remark III.3

Proposition III.2 implies that the tuning parameter $h$ in (5) can be
interpreted as the sampling time step for a forward Euler discretization of
(7) under the control input (8). Although this control interpretation is not
present in the work of Kazerounian and collaborators [9, 11, 10], it will
enable us to extend the KCM framework and obtain more realistic folding
pathways under the KCM-based simulations.

Using the KCM-based control input (8) in (7), we arrive at the following
reference vector field

$\dot{\boldsymbol{\theta}}^{\text{KCM}}:=\frac{\mathcal{J}^{\top}(\boldsymbol{\theta})\mathcal{F}(\boldsymbol{\theta})}{|\mathcal{J}^{\top}(\boldsymbol{\theta})\mathcal{F}(\boldsymbol{\theta})|_{\infty}}.$
(10)

The prominent KCM-based schemes by Kazerounian and collaborators [9, 11, 10]
for folding _in vacuo_ are therefore governed by the above reference vector
field in the protein conformation landscape, according to the forward Euler
iterations given by (5) and (6).

##  IV ODS-Based Control of Protein Folding

In this section, we will use the KCM-based reference vector field in (10) for
synthesizing QP-based control inputs for the folding dynamics in (7). The
control synthesis interpretation of protein folding enables us to add more
constraints to the folding simulation and hence making it closer to the true
folding structural biochemistry.
One important physical observation in the structural biochemistry of folding
is that to avoid extremely fast decay of excess entropy, one needs to penalize
the control input during the guided folding process (see [15, 17, 18]).
Therefore, the dihedral angle vector trajectories should not only follow the
KCM-based reference vector field in (10) but also the control torques need to
be bounded by a predetermined bound throughout the simulation. As it is shown
in “Entropy-Loss Rate of the Protein Molecule” (Section A) in Supplementary
Material, the rate of change of the protein molecule entropy linearly depends
on the control input $\mathbf{u}_{c}$ in (7). Furthermore, the upper bound on
this rate of change depends on the largest element of the control input
vector.
ODS-based control is a pointwise optimal control scheme that, in the context
of protein folding, minimizes the deviation between the closed-loop vector
field of the protein molecule dynamics in (7) and the KCM-based reference
vector field in (10), while respecting the input constraints

$|u_{c_{{}_{i}}}|\leq c_{i},\;i=1,\,\cdots,\,2N,$ (11)

where $u_{c_{{}_{i}}}$ is the $i$th control input.

In order to state the ODS-based control scheme for KCM-based protein folding,
let $\boldsymbol{\theta}(t,t_{0},\boldsymbol{\theta}_{0},\mathbf{u}_{c}(t))$,
or $\boldsymbol{\theta}(t)$ for short, denote the solution to (7) under the
control input $t\mapsto\mathbf{u}_{c}(t)$ and initial conformation
$\boldsymbol{\theta}_{0}$ at time $t_{0}$. For each solution
$\boldsymbol{\theta}(t)$, we define the set of permissible velocity vectors
$\mathcal{C}_{t}(\boldsymbol{\theta}(t))$ to be the translation by
$\boldsymbol{\theta}(t)$ of the set of free vectors

$\displaystyle\mathcal{C}(\boldsymbol{\theta}(t)):=\big{\\{}\mathbf{v}(t,\mathbf{u}_{c})\in\mathbb{R}^{2N}\big{|}\mathbf{v}=\mathcal{J}^{\top}(\boldsymbol{\theta})\mathcal{F}(\boldsymbol{\theta})+\mathbf{u}_{c},$
(12) $\displaystyle|\mathbf{u}_{c}{{}_{{}_{i}}}|\leq c_{i}\big{\\}}.$

In other words, each vector in $\mathcal{C}_{t}(\boldsymbol{\theta}(t))$
corresponds to a vector in $\mathcal{C}(\boldsymbol{\theta}(t))$ whose initial
point has been translated to $\boldsymbol{\theta}(t)$. Therefore, for any
control input $\mathbf{u}_{c}(t)$ that respects the input constraints in (11),
the possible tangent vectors to the folding pathway lie in
$\mathcal{C}_{t}(\boldsymbol{\theta}(t))$ (see Figure 2). The protein
conformation evolves according to the closed-loop vector field
$\dot{\boldsymbol{\theta}}^{\text{ODS}}$ that is not only close to the KCM-
based reference vector field $\dot{\boldsymbol{\theta}}^{\text{KCM}}$ in (10)
but also satisfies the control input constraints in (11).

Figure 2: Optimal decision strategy in the context of protein folding.

In the context of protein folding, we choose the ODS-based control law
$\mathbf{u}_{c}(t)$ in a way that at each time $t$, the instantaneous closed-
loop vector field is the “nearest” to the KCM-based reference vector field
$\dot{\boldsymbol{\theta}}^{\text{KCM}}(\boldsymbol{\theta}(t))$ given by (10)
in the norm on $\mathbb{R}^{2N}$ defined by some positive definite matrix
$\mathbf{Q}$. In other words, the control law at each $t$ is the minimizing
solution to the QP

$\displaystyle\underset{\mathbf{v}\in\mathcal{C}_{t}(\boldsymbol{\theta}(t))}{\text{min.}}\;\Big{\\{}\big{[}\mathbf{v}-\dot{\boldsymbol{\theta}}^{\text{KCM}}(\boldsymbol{\theta}(t))\big{]}^{\top}\mathbf{Q}\big{[}\mathbf{v}-\dot{\boldsymbol{\theta}}^{\text{KCM}}(\boldsymbol{\theta}(t))\big{]}\Big{\\}},$
(13)

which can be shown to be equivalent to the QP

$\displaystyle\mathbf{u}_{c}^{\ast}(\boldsymbol{\theta})=\underset{\mathbf{u}}{\text{argmin}}$
$\displaystyle\big{\\{}\frac{1}{2}\mathbf{u}_{c}^{\top}\mathbf{Q}\mathbf{u}_{c}+\mathbf{G}^{\top}(\boldsymbol{\theta})\mathbf{u}_{c}\big{\\}}$
(14) subject to $\displaystyle\;\;\;|\mathbf{u}_{c}{{}_{{}_{i}}}|\leq c_{i},$

where

$\displaystyle\mathbf{G}(\boldsymbol{\theta}):=-\mathbf{Q}\big{(}\mathcal{J}^{\top}(\boldsymbol{\theta})\mathcal{F}(\boldsymbol{\theta})-\dot{\boldsymbol{\theta}}^{\text{KCM}}(\boldsymbol{\theta})\big{)}.$
(15)

The following proposition provides sufficient conditions for uniqueness and
Lipschitz continuity of the proposed QP in (14) at each conformation
$\boldsymbol{\theta}_{0}$.

###### Proposition IV.1

Consider the QP (14), and assume the following conditions at an arbitrary
conformation $\boldsymbol{\theta}_{0}$:

  * •

The unique solution to the linear programming problem

$\displaystyle\omega=\underset{(\mathbf{u}_{c},w)\in\mathbb{R}^{2N+1}}{\text{max}}$
$\displaystyle w$ (16) subject to
$\displaystyle\;\;\;\mathbf{A}\mathbf{u}_{\text{c}}+w\mathbf{e}_{{}_{4N}}\leq\mathbf{c},$

satisfies $\omega>0$, where the vector $\mathbf{c}\in\mathbb{R}^{4N}$ is
determined by the constant bounds $c_{i}$ in (11) as
$\mathbf{c}=[c_{1},\,c_{1},\,\cdots,\,c_{2N},\,c_{2N}]^{\top}$; moreover, the
matrix $\mathbf{A}\in\mathbb{R}^{4N\times 2N}$ is given by
$\left(\begin{smallmatrix}1&0&\cdots&&0\\\ -1&0&\cdots&&0\\\
\vdots&\vdots&&&\vdots\\\ 0&&&&1\\\ 0&&&&-1\end{smallmatrix}\right),$

  * •

$\mathbf{Q}$ in (14) is symmetric and positive definite,

  * •

$\mathcal{J}^{\top}(\boldsymbol{\theta})\mathcal{F}(\boldsymbol{\theta})$ is
Lipschitz continuous at $\boldsymbol{\theta}_{0}$.

Then, the feedback $\mathbf{u}_{c}^{\ast}(\boldsymbol{\theta})$ defined in
(14) is Lipschitz continuous and unique at conformation
$\boldsymbol{\theta}_{0}$.

###### Proof:

See Section C in Supplementary Material. ∎

###### Remark IV.2

When the constraints on the control inputs in (11) are sufficiently large, or
in the limit when $c_{i}\to\infty$, the generated control inputs by the
proposed QP in (14) will be the same as the conventional KCM control torques
in (8).

##  V Simulation Results

In this section we present numerical simulation results to validate our
proposed control method. We consider a molecular chain with $N-1=10$ peptide
planes corresponding to a $2N=22$-dimensional dihedral angle space in (7). All
of our implementation has been done in MATLAB R2018b following the PROTOFOLD I
guidelines [9, 11] on an Intel® Core™ i7-6770HQ<EMAIL_ADDRESS>In our current
implementation, the side chains are ignored. See Section E of Supplementary
Material for the implementation details and computational complexity analysis
of the overall algorithm.

We run three sets of simulations. Our simulations are based on the forward
Euler discretization of the underlying continuous time dynamics in (7).
Propositions III.2 and IV.1 guarantee the “well-behavedness” of the resulting
trajectories under small enough time steps and not too tight constraints on
the control inputs. The first simulation is run by using the original KCM-
based control torque in (8) with $h=0.04$ and for 325 iterations starting from
a vector of initial dihedral angles uniformly distributed about the mean value
$27.7^{{}^{\circ}}$ with standard deviation $1.1^{{}^{\circ}}$111See
Supplementary Material for the exact value of $\boldsymbol{\theta}_{0}$. . The
time it takes to run the program is equal to 3.47 sec. As expected from the
KCM-based folding scheme [9, 10], the torques converge to a small neighborhood
of zero and the free energy of the molecule also converge to a local minimum
on the folding landscape. The blue curves in 3(a) and 3(b) depict the profiles
of the maximum control torque $|\mathbf{u}_{c}|_{\infty}$ and the aggregated
free energy $\mathcal{G}(\boldsymbol{\theta})$ of the conventional KCM
framework versus the number of iterations, respectively. Next, we constrain
the control torques used for guiding the folding process according to the QP-
based control scheme in (14) while using the KCM reference vector field in
(10). We use

$\displaystyle\mathbf{Q}=\sqrt{2N}\mathbf{I}_{2N},\;c_{i}=c_{0}\frac{\rho}{\sqrt{2N}},\,1\leq
i\leq 2N,$ (17)

for the quadratic weight matrix and the constraints on the folding control
torques, respectively. In (17), the constant factor $c_{0}$ for conversion to
$\text{kcal}/\text{mol}$ units of energy is equal to
$N_{0}q_{e}^{2}k_{\text{J-kcal}}\times 10^{20}$ where $N_{0}=6.022\times
10^{23}$ is the Avogadro number, $q_{e}=1.602\times 10^{-19}$ is the charge of
electron in coulombs, $k_{\text{J-kcal}}=1/4184$ is the energy unit conversion
factor from Joules to kilo calories. The parameter $\rho$ in (17) is used for
relaxing/tightening the constraints on the control torques.

(a) Joint torque profiles during protein folding.

(b) Energy profiles during protein folding.

(c) Initial and final conformations.

Figure 3: Simulation results using the conventional and the QP-based KCM
control torques with two different sets of constraints.

The orange (associated with $\rho=20$) and red (associated with $\rho=9$)
curves in 3(a) and 3(b) depict the profiles of $|\mathbf{u}_{c}|_{\infty}$ and
the free energy of the peptide chain $\mathcal{G}(\boldsymbol{\theta})$ versus
the number of iterations, respectively. As it can be seen from 3(a), the QP-
based control torques respect the predetermined bounds while both converging
to zero. It takes around 4.1sec to run each of these two simulations. This
should not be surprising as we are running a QP to respect the control input
constraints during folding (see Section E of Supplementary Material). Indeed,
the increased computational complexity is the price to be paid for obtaining
more realistic folding pathways. Next, to violate the first condition in
Proposition 4.1, we gradually made the constraints tighter to the point of
failure at around $\rho=2$, where QP-based control fails to exist at certain
conformations along the folding pathway. While both KCM and QP-based KCM
converge to the same local minimum of the protein molecule energy, the folding
pathways generated by the QP-based KCM result in both smaller control torques
and transient energy values; hence, respecting the biological observations
about the entropy-loss constraints during folding.

##  VI Concluding Remarks and Future Research Directions

Considering the importance of computer-aided prediction of the folded
structure of protein molecules and their folding pathways, this paper formally
extended the KCM-based framework for protein conformation prediction by
including entropy-loss constraints via a QP-based control solution. Indeed,
constraining the entropy-loss rates during the folding simulations would make
the obtained pathways more realistic. After relating the conventional KCM
iterations to a proper nonlinear affine control system, we demonstrated that
the well-established KCM torques can be utilized for defining a reference
vector field for QP-based control schemes that account for the entropy-loss
constraints on the folding pathway. Furthermore, we provided sufficient
conditions for uniqueness and Lipschitz continuity of the generated control
torques at each conformation of the protein molecule. Simulations for a
simplified molecular chain demonstrated the effectiveness of the proposed
methodology in obtaining more realistic folding pathways. These findings lead
us to further research avenues such as investigating the folding pathways for
two-state proteins, protein misfolding dynamics, and developing adaptive
schemes for TMD-based folding.

## References

  * [1] G. S. Chirikjian, K. Kazerounian, and C. Mavroidis, “Analysis and design of protein based nanodevices: Challenges and opportunities in mechanical design,” _J. Mech. Des._ , vol. 127, no. 4, pp. 695–698, 2005.
  * [2] A. V. Finkelstein and O. Ptitsyn, _Protein physics: a course of lectures_. Elsevier, 2016.
  * [3] F. Bergasa-Caceres and H. A. Rabitz, “Interdiction of protein folding for therapeutic drug development in SARS CoV-2,” _J. Phys. Chem. B_ , vol. 124, no. 38, pp. 8201–8208, 2020.
  * [4] M. AlQuraishi, “A watershed moment for protein structure prediction,” _Nature_ , vol. 577, pp. 627–628, 2020.
  * [5] J. Moult, K. Fidelis, A. Kryshtafovych, T. Schwede, and A. Tramontano, “Critical assessment of methods of protein structure prediction (CASP)–round x,” _Proteins: Struct. Funct. Bioinform._ , vol. 82, pp. 1–6, 2014.
  * [6] S. Ołdziej, C. Czaplewski, A. Liwo, M. Chinchio, M. Nanias, J. Vila, M. Khalili, Y. Arnautova, A. Jagielska, M. o. Makowski _et al._ , “Physics-based protein-structure prediction using a hierarchical protocol based on the UNRES force field: assessment in two blind tests,” _Proc. Natl. Acad. Sci._ , vol. 102, no. 21, pp. 7547–7552, 2005.
  * [7] R. A. Lippert, K. J. Bowers, R. O. Dror, M. P. Eastwood, B. A. Gregersen, J. L. Klepeis, I. Kolossvary, and D. E. Shaw, “A common, avoidable source of error in molecular dynamics integrators,” _J. Chem. Phys._ , vol. 126, no. 4, p. 046101, 2007.
  * [8] K. Kazerounian, K. Latif, K. Rodriguez, and C. Alvarado, “Nano-kinematics for analysis of protein molecules,” _J. Mech. Des._ , vol. 127, no. 4, pp. 699–711, 2005.
  * [9] K. Kazerounian, K. Latif, and C. Alvarado, “Protofold: A successive kinetostatic compliance method for protein conformation prediction,” _J. Mech. Des._ , vol. 127, no. 4, pp. 712–717, 2005.
  * [10] P. Tavousi, M. Behandish, H. T. Ilieş, and K. Kazerounian, “Protofold II: Enhanced model and implementation for kinetostatic protein folding,” _J. Nanotechnol. Eng. Med._ , vol. 6, no. 3, 2015.
  * [11] P. Tavousi, “On the systematic design and analysis of artificial molecular machines,” Ph.D. dissertation, University of Connecticut, 2016.
  * [12] Z. Shahbazi, H. T. Ilieş, and K. Kazerounian, “Hydrogen bonds and kinematic mobility of protein molecules,” _J. Mech. Robot._ , vol. 2, no. 2, 2010.
  * [13] C. Madden, P. Bohnenkamp, K. Kazerounian, and H. T. Ilieş, “Residue level three-dimensional workspace maps for conformational trajectory planning of proteins,” _Int. J. Robot. Res._ , vol. 28, no. 4, pp. 450–463, 2009\.
  * [14] D. M. Dunlavy, D. P. O’leary, D. Klimov, and D. Thirumalai, “HOPE: A homotopy optimization method for protein structure prediction,” _J. Comput. Biol._ , vol. 12, no. 10, pp. 1275–1288, 2005.
  * [15] Y. Arkun and B. Erman, “Prediction of optimal folding routes of proteins that satisfy the principle of lowest entropy loss: dynamic contact maps and optimal control,” _PloS One_ , vol. 5, no. 10, p. e13275, 2010.
  * [16] Y. Arkun and M. Gür, “Combining optimal control theory and molecular dynamics for protein folding,” _PLoS One_ , vol. 7, no. 1, p. e29628, 2012\.
  * [17] ——, “Protein folding using coarse-grained optimal control and molecular dynamics,” _IFAC Proceedings Volumes_ , vol. 44, no. 1, pp. 14 213–14 216, 2011.
  * [18] T. R. Weikl and K. A. Dill, “Folding rates and low-entropy-loss routes of two-state proteins,” _J. Mol. Biol._ , vol. 329, no. 3, pp. 585–598, 2003\.
  * [19] J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kale, and K. Schulten, “Scalable molecular dynamics with NAMD,” _J. Comput. Chem._ , vol. 26, no. 16, pp. 1781–1802, 2005.
  * [20] R. Barnard, “An optimal-aim control strategy for nonlinear regulation systems,” _IEEE Trans. Automat. Contr._ , vol. 20, no. 2, pp. 200–208, 1975\.
  * [21] M. Spong, J. Thorp, and J. Kleinwaks, “The control of robot manipulators with bounded input,” _IEEE Trans. Automat. Contr._ , vol. 31, no. 6, pp. 483–490, 1986.
  * [22] A. Mohammadi and M. W. Spong, “Integral line-of-sight path following control of magnetic helical microswimmers subject to step-out frequencies,” _Automatica_ , vol. 128, p. 109554, 2021.
  * [23] R. Thomas, J. Thorp, and C. Pottle, “A model-referenced controller for stabilizing large transient swings in power systems,” _IEEE Trans. Automat. Contr._ , vol. 21, no. 5, pp. 746–750, 1976.
  * [24] V. Murali, N.-S. P. Hyun, and E. I. Verriest, “Graceful transitions between periodic walking gaits of fully actuated bipedal robots,” in _2020 American Contr. Conf. (ACC)_ , 2020, pp. 1043–1048.
  * [25] V. Klemm, A. Morra, L. Gulich, D. Mannhart, D. Rohr, M. Kamel, Y. de Viragh, and R. Siegwart, “LQR-assisted whole-body control of a wheeled bipedal robot with kinematic loops,” _IEEE Robot. Autom. Lett._ , vol. 5, no. 2, pp. 3745–3752, 2020.
  * [26] A. D. Ames, X. Xu, J. W. Grizzle, and P. Tabuada, “Control barrier function based quadratic programs for safety critical systems,” _IEEE Trans. Automat. Contr._ , vol. 62, no. 8, pp. 3861–3876, 2017.
  * [27] A. G. Pandala, Y. Ding, and H.-W. Park, “qpswift: A real-time sparse quadratic program solver for robotic applications,” _IEEE Robot. Autom. Lett._ , vol. 4, no. 4, pp. 3355–3362, 2019.
  * [28] A. Amadei, A. B. Linssen, and H. J. Berendsen, “Essential dynamics of proteins,” _Proteins: Struct. Funct. Bioinform._ , vol. 17, no. 4, pp. 412–425, 1993.
  * [29] Q. Shao, H. Wei, and Y. Q. Gao, “Effects of turn stability and side-chain hydrophobicity on the folding of $\beta$-structures,” _J. Mol. Biol._ , vol. 402, no. 3, pp. 595–609, 2010.
  * [30] C. Alvarado and K. Kazerounian, “On the rotational operators in protein structure simulations,” _Prot. Eng._ , vol. 16, no. 10, pp. 717–720, 2003\.
  * [31] D. B. Adolf and M. D. Ediger, “Brownian dynamics simulations of local motions in polyisoprene,” _Macromolecules_ , vol. 24, no. 21, pp. 5834–5842, 1991\.
  * [32] P. Ferrara, J. Apostolakis, and A. Caflisch, “Computer simulations of protein folding by targeted molecular dynamics,” _Proteins: Struct. Funct. Bioinform._ , vol. 39, no. 3, pp. 252–260, 2000.

## Supplementary Material To “Quadratic Optimization-Based Nonlinear Control
for Protein Conformation Prediction”

Notation. In these supplementary notes, we will number the equations using
(s#), the references using [S#], and the figures using SF#. Therefore, (s1)
refers to Equation 1 in the supplementary notes while (1) refers to Equation 1
in the original article. The same holds true for the references and the
figures.

###  VI-A Entropy-Loss Rate of the Protein Molecule

There are a few equivalent ways to compute the entropy changes of a protein
molecule during folding including the methods based on configurational entropy
computations and the methods based on thermodynamics-based operational
definitions. Here, we will use Clausius thermodynamics expression in [S1],
which is similar to the derivations in [15] in that it is based on classical
thermodynamics arguments. Based on Clausius equation, the change of the
molecule entropy is given by

$\Delta
S(\boldsymbol{\theta})=\kappa_{0}\Delta\mathcal{G}(\boldsymbol{\theta}),$ (s1)

where $\kappa_{0}$ is a temperature-dependent constant and
$\mathcal{G}(\cdot)$ is the aggregated free energy of the molecule given by
(2). Hence,

$\frac{d}{dt}S(\boldsymbol{\theta})=\kappa_{0}\frac{\partial\mathcal{G}}{\partial\boldsymbol{\theta}}\dot{\boldsymbol{\theta}}=\kappa_{0}\frac{\partial\mathcal{G}}{\partial\boldsymbol{\theta}}\big{(}\mathcal{J}^{\top}(\boldsymbol{\theta})\mathcal{F}(\boldsymbol{\theta})+\mathbf{u}_{c}\big{)}.$
(s2)

Using Hölder’s and the triangle inequalities in (s2), it can be seen that

$|\dot{S}|\leq\kappa_{0}\bigg{|}\frac{\partial\mathcal{G}}{\partial\boldsymbol{\theta}}\bigg{|}_{1}\bigg{(}\big{|}\mathcal{J}^{\top}(\boldsymbol{\theta})\mathcal{F}(\boldsymbol{\theta})\big{|}_{\infty}+\big{|}\mathbf{u}_{c}\big{|}_{\infty}\bigg{)}.$
(s3)

where $|\cdot|_{1}$ denotes the vector 1-norm.

###  VI-B In-Depth Comparison with the Optimal Control Framework in [15, 17,
16]

One of the inspirations for our ODS-based control scheme for protein folding
in (14) has been due to the pioneering work by Arkun and collaborators [15,
17, 16]. Here, we provide an in-depth comparison with our work and state
possible extensions of the work by Arkun and collaborators [15, 17, 16] to our
nonlinear setting.

  * •

Difference in the studied dynamics: As we stated in Remark 3.1 of the article,
we are considering a _nonlinear dynamics_ (as opposed to the _linear dynamics_
in [15, 17, 16]) in the dihedral angle space of the molecule (as opposed to
[15, 17, 16], where the state vector of the system is the Cartesian position
of the $\text{C}_{\alpha}$ atoms).

  * •

Nonlinear extension of [15, 17, 16]: In the linear time-invariant (LTI) setup
considered by Arkun and collaborators [15, 17, 16], the control input is
penalized by using an infinite horizon LQR cost function to avoid extremely
fast decay of excess entropy. In the context of the KCM framework, a cost
function in the dihedral angle space can be written similar to [15, 17, 16] as

$J[\mathbf{u}_{c}]:=\int_{0}^{T_{f}}\big{(}\tilde{\boldsymbol{\theta}}(t)^{\top}\mathbf{Q}\tilde{\boldsymbol{\theta}}(t)+\rho\mathbf{u}_{c}(t)^{\top}\mathbf{P}\mathbf{u}_{c}(t)\big{)}dt,$
(s4)

where
$\tilde{\boldsymbol{\theta}}:=\boldsymbol{\theta}-\boldsymbol{\theta}^{\text{ref}}$
is the error from the current to a reference conformation, e.g., the native
conformation of the protein. Furthermore, the tunable parameter $\rho>0$ is
employed to make a trade-off between avoiding the high energy regions of the
landscape while choosing entropically preferred folding pathways by penalizing
high-entropy-loss routes. The similarities with [15, 17, 16] stop here. Under
an LTI control setup (which is not the case for the nonlinear dynamics in (7))
and letting $T_{f}\to\infty$, the solution to (s4) takes the form of the
feedback control law
$\mathbf{u}_{c}^{\ast}(\boldsymbol{\theta})=-\mathbf{K}(\rho)\tilde{\boldsymbol{\theta}}$,
where the state feedback gain matrix $\mathbf{K}(\rho)$ is the solution to a
proper algebraic Riccati equation (ARE). However, this ARE-based feedback
control law is not the minimzer of the cost function in (s4) when the dynamics
are governed by (7).
Indeed, if the system has nonlinearities, which is the case in this article,
there is a need for solving a Hamilton-Jacobi equation (HJE). HJEs are partial
differential equations whose explicit closed-form solutions do not exist, even
in the case of a simple nonlinearity (see, e.g., [S3]).
Hamilton-Jacobi equation (HJE) for the protein folding dynamics given by (7)
can obtained as follows. Let us consider the nonlinear dynamics in (7) and the
cost function in (s4). Defining
$\tilde{\boldsymbol{\theta}}:=\boldsymbol{\theta}-\boldsymbol{\theta}^{\text{ref}}$,
$\tilde{\mathbf{G}}(\tilde{\boldsymbol{\theta}}):=\mathcal{J}^{\top}(\tilde{\boldsymbol{\theta}}+\boldsymbol{\theta}^{\text{ref}})\mathcal{F}(\tilde{\boldsymbol{\theta}}+\boldsymbol{\theta}^{\text{ref}})$,
the optimal control input is given by [S4]

$\mathbf{u}_{c}^{\ast}(\tilde{\boldsymbol{\theta}})=-\frac{1}{2\rho}\mathbf{P}^{-1}\nabla
V(\tilde{\boldsymbol{\theta}})+\tilde{\mathbf{G}}(\mathbf{0}).$ (s5)

where $V(\cdot)$ is a twice continuously differentiable function that
satisfies HJE

$\frac{1}{4\rho}\nabla V^{\top}\mathbf{P}^{-1}\nabla
V-(\mathcal{J}^{\top}(\boldsymbol{\theta})\mathcal{F}(\boldsymbol{\theta})-\tilde{\mathbf{G}}(\mathbf{0}))\nabla
V-\tilde{\boldsymbol{\theta}}^{\top}\mathbf{Q}\tilde{\boldsymbol{\theta}}=0.$
(s6)

For solving the above numerically demanding partial differential equation,
there are a few methods such as employing artificial neural networks (see,
e.g., [S4]).

  * •

Computational load: Even if we neglect the complexity arising from solving
HJEs in a nonlinear setting, the framework in [15, 17, 16] requires _at least
twice the computational time_ of our framework. This is due to the fact that
the framework in [15, 17, 16] requires knowing the final folded structure of
the protein to find the proper folding pathways. Our work, on the other hand,
does not require knowing this final folded structure. Rather, by following the
reference vector field coming from the KCM framework, folding starts from a
random structure and converges to a local minimum of the aggregate free energy
of the molecule. Therefore, a fair comparison with the nonlinear extension of
[15, 17, 16] would be to run the simulations in the dihedral angle space
twice; namely, finding the final folded structure
$\boldsymbol{\theta}^{\text{ref}}$ using the conventional KCM simulation in
the first run, and, in the second run, using the nonlinear extension of the
optimal control framework in [15, 17, 16] (as formulated in (s5), (s6)) to
guide the molecule towards the final conformation
$\boldsymbol{\theta}^{\text{ref}}$, which has been obtained from the first
run.

  * •

The underlying philosophy of the pointwise optimal control laws: ODS-based
control and its resulting QP-based synthesis [21] were originally developed to
address some of the challenges arising from integral performance measures such
as the one in (s4). Indeed, the QP-based synthesis avoids problems such as
computationally intensive schemes and bypassing the Pontryagin maximum
principle and/or HJEs under the bound constraints on the allowable input
torques. Even ignoring the computational complexity of such schemes, one might
prefer an instantaneous or pointwise optimization procedure which optimizes
the present state of the system without regard to future events. In the
context of our QP-based control framework, the reference vector field in (10)
guides the protein towards its native conformation on the folding landscape.

###  VI-C Proof of Proposition IV.1

###### Proof:

The solution $\omega$ to the linear programming problem in (16) is the width
of a feasible set associated with the QP in (14). Furthermore, Lipschitz
continuity of
$\mathcal{J}^{\top}(\boldsymbol{\theta})\mathcal{F}(\boldsymbol{\theta})$ at
conformation $\boldsymbol{\theta}_{0}$ implies Lipschitz continuity of
$\mathbf{G}(\boldsymbol{\theta})$ in the QP. Following Theorem 1 in [S5],
which is based on Mangasarian-Fromovitz regularity conditions [S6], uniqueness
and Lipschitz continuity of the control input follow. ∎

###  VI-D Proof of Proposition III.2

###### Proof:

Employing the control input (8) in (7) and noting that the iteration in (5) is
the forward Euler difference equation for the closed-loop dynamics with time
step $h$, the rest of the proof follows a standard argument from numerical
analysis (see, e.g., the proof of Theorem 1.1 in [S7]). ∎

###  VI-E Details of Implementation and Computational Complexity of the
Overall Algorithm

Details of implementation. A flowchart of the program is depicted in Figure
SF1. The convergence criteria, following PROTOFOLD I [9], is that the
magnitude of all of the joint equivalent torques to be less than or equal to 1
kcal/(mol.A), which is in accordance with the experimentally determined values
in folded proteins. Our QP-based control input generation block is using the
MATLAB ‘quadprog’ routine, where we are using the ‘interior-point-convex’ as
the QP numerical solver. After performing the direct kinematics calculations
on the chain, we are using the ‘molecule3D’ program [S9] for drawing the
protein molecules. Interested readers are encouraged to contact the
corresponding author<EMAIL_ADDRESS>for obtaining a copy of the
program.

Computational complexity. The big O notation for complexity analysis of
numerical algorithms can be defined as follows. Given two real-valued
functions $f(\cdot)$ and $g(\cdot)$, we say that $f(x)=\mathcal{O}(g(x))$ if
there exists a positive real number $M$ and a real number $x_{0}$ such that
$|f(x)|\leq M|g(x)|$ for all $x\geq x_{0}$.

The QP given by (14) belongs to the class of box-constrained convex quadratic
programs that can be efficiently solved using interior-point algorithms (see,
e.g., [S8]). The time needed for solving such a QP using interior-point
methods is of order $\mathcal{O}(N^{3})$, where there are $N-1$ peptide planes
in the protein molecule. Furthermore, if there are $N_{a}$ atoms in the chain,
the computational complexity associated with computing the intramolecular
forces at each iteration is of order $\mathcal{O}(N_{a}^{2})$. Finally, the
computational complexity associated with direct kinematics calculations is of
order $\mathcal{O}(N)$.

Figure SF1: The flowchart of the computer program. Contact the corresponding
author<EMAIL_ADDRESS>for obtaining a copy of the program.

Initial conformation of the protein molecule. In the presented simulation
results in the paper, the initial protein conformation vector (in degrees) is
given by

$\displaystyle\boldsymbol{\theta}_{0}=\big{[}$ $\displaystyle
28.3,28.6,26.1,28.7,27.7,26.0,26.6,27.5,28.8,\cdots$ $\displaystyle
28.8,26.2,28.8,28.8,27.2,28.3,26.2,27.1,28.7,\cdots$ $\displaystyle
28.3,28.8,27.8,25.8\big{]}^{\circ\top}.$

Figure SF2: The protein molecule kinematic structure under the
$\text{C}_{\alpha}-\text{CO}-\text{NH}-\text{C}_{\alpha}$ coplanarity
assumption.

## Supplementary References

[S1] S. Hikiri, T. Yoshidome, and M. Ikeguchi. “Computational methods for
configurational entropy using internal and Cartesian coordinates,” _J. Chem.
Theory Comput._ , vol. 12, no. 12, pp. 5990-6000, 2016.

[S2] T. Ohtsuka. “Solutions to the Hamilton-Jacobi equation with algebraic
gradients,” _IEEE Trans. Automat. Contr._ , 56(8), pp. 1874-1885, 2010.

[S3] C-H. Won, S. Biswas. “Optimal control using an algebraic method for
control-affine non-linear systems,” Int. J. Contr., vol. 80, no. 9, pp.
1491-1502, 2007.

[S4] T. Cheng, F. L. Lewis and M. Abu-Khalaf. “Fixed-final-time-constrained
optimal control of nonlinear systems using neural network HJB approach,” _IEEE
Trans. Neural Netw._ , vol. 18, no. 6, pp.1725-1737, 2007.

[S5] B. Morris, M. J. Powell, and A. D. Ames. “Sufficient conditions for the
Lipschitz continuity of QP-based multi-objective control of humanoid robots,”
in _Proc. 52nd IEEE Conf. Dec. Contr._ , Firenze, Italy, Dec. 2013, pp.
2920-2926.

[S6] O. L., Mangasarian, and S. Fromovitz. “The Fritz John necessary
optimality conditions in the presence of equality and inequality constraints,”
J. Math. Anal. Appl., vol. 17, no. 1, pp. 37-47, 1967.

[S7] A. Iserles, _A first course in the numerical analysis of differential
equations_. Cambridge University Press, 2009, no. 44.

[S8] Y. Wang, and S. Boyd. “Fast evaluation of quadratic control-Lyapunov
policy,” IEEE Trans. Contr. Syst. Technol., vol. 19, no. 4, pp. 939-946, 2010.

[S9] A. Ludwig. _molecule3D_ (2016). Accessed: Oct. 2020. [Online]. Available:
https://www.mathworks.com/matlabcentral/fileexchange/55231-molecule3d, MATLAB
Central File Exchange.