# Single-grasp deformable object discrimination: the effect of gripper
morphology, sensing modalities, and action parameters

Michal Pliska1, Michal Mares1, Pavel Stoudek1, Zdenek Straka1, Karla
Stepanova2, and Matej Hoffmann1 1Department of Cybernetics, Faculty of
Electrical Engineering, Czech Technical University in Prague;
<EMAIL_ADDRESS>. 2Czech Institute of Informatics, Robotics, and
Cybernetics, Czech Technical University in Prague; This work was supported by
the project Interactive Perception-Action-Learning for Modelling Objects
(IPALM) (H2020 – FET – ERA-NET Cofund – CHIST-ERA III / Technology Agency of
the Czech Republic, EPSILON, nr. TH05020001). M.H. and Z.S. were also
partially supported by OP VVV MEYS funded project
CZ.02.1.01/0.0/0.0/16_019/0000765 “Research Center for Informatics”. K.S. was
partially supported by the Czech Science Foundation (project no. GA21-31000S).

###### Abstract

We studied the discrimination of deformable objects by grasping them using 4
different robot hands / grippers: Barrett hand (3 fingers with adjustable
configuration, 96 tactile, 8 position, 3 torque sensors), qb SoftHand (5
fingers, 1 motor, position and current feedback), and two industrial type
parallel jaw grippers with position and effort feedback (Robotiq 2F-85 and
OnRobot RG6). A set of 9 ordinary objects differing in size and stiffness and
another highly challenging set of 20 polyurethane foams differing in material
properties only was used. We systematically compare the grippers’ performance,
together with the effects of: (1) type of classifier (k-NN, SVM, LSTM)
operating on raw time series or on features, (2) action parameters (grasping
configuration and speed of squeezing), (3) contribution of sensory modalities.
Classification results are complemented by visualization of the data using
PCA. We found: (i) all the grippers but the qb SoftHand could reliably
distinguish the ordinary objects set; (ii) Barrett Hand reached around 95%
accuracy on the foams; OnRobot RG6 around 75% and Robotiq 2F-85 around 70%;
(iii) across all grippers, SVM over features and LSTM on raw time series
performed best; (iv) faster compression speeds degrade classification
performance; (v) transfer learning between compression speeds worked well for
the Barrett Hand only; transfer between grasping configurations is limited;
(vi) ablation experiments provided intriguing insights—sometimes a single
sensory channel suffices for discrimination. Overall, the Barrett Hand as a
complex and expensive device with rich sensory feedback provided best results,
but uncalibrated parallel jaw grippers without tactile sensors can have
sufficient performance for single-grasp object discrimination based on
position and effort data only. Transfer learning between the different robot
hands remains a challenge.

##  I Introduction

Although visual object recognition has seen tremendous progress in the last
years, not all object properties can be perceived using distal sensing. In
particular, physical object properties like stiffness, roughness, or mass are
better acquired through manipulation. Object recognition from haptic
exploration can be also more robust as it is insensitive to lighting
conditions and object attributes that may be irrelevant for a task (e.g.,
color). Deformable objects constitute an important case. Their softness or
elasticity can hardly be perceived without making contact with the object.
Furthermore, gauging their specific material properties may be the goal such
as when estimating fruit ripeness.

Sanchez et al. [1] provide a survey of robotic manipulation and sensing of
deformable objects. In this work, we focus on triparametric objects—solid
objects such as sponges or plush toys, which are also the least researched
object type [1]. Our motivation is practical. Unlike works that employ robot
hands with tactile arrays and complex haptic exploratory procedures (e.g., [2,
3]), we focus on single-grasp recognition—a “functional action” [4] that
serves an application goal. The object is picked up by a manipulator and
grabbing of the object itself should immediately serve to classify the object
such that it can be sorted like fruit based on its ripeness or
plastic/paper/metal on a single-stream recycling line [5, 6].

Computational costs when simulating interactions with deformable objects such
as using the Finite Element Method (FEM) are enormous and prevent online
recognition (see [7] for a review). The alternative are data-driven methods
that directly process the time series generated by the interaction with the
object. In this case, the detailed characteristics of the sensory streams
become critical. This involves not only their number and type (force,
effort/current, touch, encoders/position sensors, etc.) but also the way they
are generated—the gripper morphology and the parameters of the grasping action
profoundly shape what is perceived (see morphology facilitating perception in
[8] for a theoretical account and [9] for applying this idea to terrain
recognition by a running robot).

Contribution. First, two sets of deformable objects were grasped using 4
different robot hands / grippers: Barrett hand (3 fingers with adjustable
configuration), qb SoftHand (5 fingers, 1 motor), and two industrial type
parallel jaw grippers (Robotiq 2F-85 and OnRobot RG6). To our knowledge, other
researchers have employed a single robot hand or gripper. The dataset is made
publicly available. Second, depending on the gripper type, we analyze the
effect of grasp parameters—grasping configuration and compression speed—and
the performance of transfer learning. Third, we study the effect of individual
sensory channels by performing object classification with a subset of the
sensory channels (ablation). Fourth, we trained 4 different classifiers
(feature-based k-NN and SVM; feature-less k-NN and LSTM) and demonstrate that
that most of the observations regarding the effect of gripper morphology,
sensory channels, and action parameters generalize across the input type and
classifier. Fifth, object classification is accompanied by unsupervised
analysis of the feature space using dimensionality reduction techniques (PCA),
allowing to assess the relative effects of the parameters under investigation:
gripper, action parameters, and object. An interactive tool for additional
visualizations is publicly available. Finally, the dataset of 20 polyurethane
foams is to our knowledge more challenging than other datasets that have been
studied so far. Nevertheless, we show that single-grasp discrimination using
uncalibrated sensors with a relatively low sampling rate is successful with 3
of the 4 devices tested.

##  II Related work

Recent surveys of haptic or tactile robot perception are provided by [10, 11].
For material recognition by tactile sensing, Luo et al. [10] distinguish
surface texture and object stiffness based recognition. The former, not
addressed in this work, typically involves sliding along the object surface
and perceiving the vibrations using acoustic or tactile sensors. The latter,
object stiffness based tactile material recognition, may involve tapping on
the object, pressing it against a surface, or squeezing it between the gripper
jaws. Tapping is mainly estimating the object’s hardness, while pressing or
squeezing probes its stiffness or material elasticity. If multiple modalities
can be perceived such as when using the BioTac sensors (normal and shear
force, vibrations, and temperature) [3, 12, 2] or the skin modules measuring
normal force, proximity sensor, acceleration, and temperature [13], it becomes
possible to discriminate objects based on a combination of these properties.
Li et al. [11] distinguish contact-level information from individual contact
sites and object-level information regarding objects as a whole, often
acquired from multiple contacts.

The focus of this work is not tactile object recognition as an endeavor
exploiting several submodalities of touch (normal and tangential force,
vibration, temperature) and involving exploration of the object’s surface.
Instead, we focus on compressing the object between the fingers of a parallel
jaw gripper or a robot hand and using only the mechanical response from
squeezing the object—the stress/strain curve—for discrimination. Only one of
the devices we used, the Barrett Hand, was equipped with tactile sensors and
thus we can assess their specific contribution as well.

Two-finger grippers are the simplest devices. Chin et al. [5] used a custom-
made soft two-finger gripper with pressure and strain sensors. An
underactuated compliant two-finger gripper with 8 tactile sensors on every
finger was used in [14, 4]. More standard parallel jaw grippers with tactile
sensors were used in [6, 15, 16].

Multi-finger hands allow for more sophisticated manipulation, enabling in-hand
exploratory sequences. Delgado et al. [17] used the Shadow hand to grasp
deformable objects and computed the deformability ratio to infer the maximum
force allowed to be exerted on the object. Object classification was performed
in [18, 19]. A 3-finger hand and a parallel jaw gripper were used in [20].

###  II-A Sensory channels

A line can be drawn between devices that are capable of multimodal sensing
(e.g., [3, 12, 2, 13]), tactile arrays which are composed of several sensors
perceiving normal force/pressure [21, 22, 6, 14, 15, 4], and perceiving the
force-displacement (stress-strain) relationship only through the gripper
(typically from the motor effort). Tactile sensing from fingertips may be
combined with motor positions and currents [19]. For parallel jaw grippers,
one tactile array may be put on every gripper jaw [15]; sometimes a single
pressure sensor is attached to every jaw [16]; Guo et al. [6] used a tactile
array on only one gripper jaw.

###  II-B Features

We primarily focus on object compression between the jaws of a parallel
gripper, but other haptic exploration actions like squeezing an object between
fingers with tactile sensors share some characteristics. Typically, the
compression action gives rise to a time series with a minimum of two channels:
a position channel (like distance between gripper jaws) and force channel
(resistance perceived by the gripper or pressure sensed by tactile sensors).
There may be more channels. For position, in case a hand with fingers is used,
joint angles from encoders may be available. For force, there may be many
readings from individual tactile/pressure sensors.

The most common approach in the literature has been to compute features from
these time series. There are two distinct events: (i) the moment of contact
with the object and (ii) the final position when compression is stopped. In
[14, 4], these events were explicitly exploited, using sensor readings only at
these instances. Others used features calculated over the whole time series.
The features may be physically motivated like the change in force/pressure
over the compression (change in length) [2] or max. force, contact area, and
contact motion in [21] or more generic (mean, variance, maximum, minimum,
range features in [6]). For multimodal data like from the BioTac sensors, more
complex feature sets can be constructed [3, 12]. Dimensionality reduction
techniques like Principal Component Analysis (PCA) can be applied on top of
the features [21] or on raw time series data [12]. Specific features to encode
the time series from tactile sensors (Spatio-Temporal Hierarchical Matching
Pursuit [20]) or force-displacement curve from a gripper (functional PCA [16])
can be designed.

Recurrent neural networks can be fed directly with time series data, as used
in [18] to feed a Jordan-type network. This approach may gain in popularity
with the recent advent of classifiers based on LSTMs.

###  II-C Classifiers

Several classifiers were employed to classify objects based on deformation
cues. Typically, features were used as inputs. The most popular classifiers
were Support Vector Machines (SVM) [6, 12, 23, 14, 16], k-nearest neighbor
(k-NN) [21, 22, 16], random forests [12, 14], naïve Bayes [12, 16], Hidden
Markov Models [22], feed-forward neural networks like Multi-Layer Perceptron
[19, 16], and recurrent neural networks—Jordan type [18] or, more recently,
LSTMs [24, 22, 19].

###  II-D Exploratory actions

Sometimes, humans rely on only brief haptic exposures, what Klatzky and
Lederman called a haptic glance [25]. In robotics, an analogue of such a
haptic glance may be a single grasp [4]. The action parameters of such a
grasp—its speed and maximum force—are rarely studied. Often, open loop control
with a given velocity is applied until a stopping condition, which can be a
force threshold or that the position has reached a steady state. Similarly to
the concepts of stiffness and elasticity in physics, compression is the
independent variable and the resistance to deformation/force/pressure the
dependent variable. The maximum force threshold should be set correctly in
order to induce sufficient compression while at the same time not damaging the
object (compressing beyond the elastic limit)—investigated in [16]. However,
some works take the opposite approach maintaining a constant force and
measuring deformation [2]. Different compression velocities were tested in
[16]. Bhattacharjee et al. [22], using skin on a robot forearm, studied the
effect of velocity and stiffness of the manipulator as a control parameter on
classification of the objects the arm collided with.

##  III Experimental setup

###  III-A Objects set

We considered two sets of deformable objects. The first set, ordinary objects
set, consisted of 9 mostly cuboid objects with a different size and degree of
deformability. The second set, polyurethane foams set, consisted of 20
polyurethane foam blocks of similar size, along with reference values for
elasticity and density provided by the manufacturer. Cuboids are preferred
over spheres, as the contact surface area does not change during deformation.
For all the objects used, the deformation can be regarded as elastic—objects
returning to their undeformed shapes once the external force is
removed—although for some objects, this recovery was slow (e.g., Blue die and
Blue cube).

####  III-A1 Ordinary objects set

This set of 9 deformable objects is visualized in Fig. 1(a), approximately
spread out by the stiffness/elasticity of the objects and their volume. The
_yellowcube_ is composed of the same material as the _yellowsponge_ —it has
been cut out from another exemplar of the same sponge, aiming at the dimension
of the Kinova cube. The same is true for _bluedie_ and _bluecube_. Conversely,
_whitedie_ , _kinova cube_ , _yellowcube_ and _bluecube_ have roughly same
dimensions but different material composition and hence stiffness. The dataset
has been deliberately designed in this way, in order to test which of the
object properties are key for model-free haptic object recognition.

(a)

(b)

Figure 1: (a) Ordinary objects set approximately spread out on the elasticity
and volume axes (reference values for this object set are not available). (b)
Polyurethane foams set.

####  III-A2 Polyurethane foams set

To complement the “ordinary objects set”, we tested on another set of objects:
20 polyurethane foams (Fig. 1(b)). These were samples provided by a
manufacturer of mattresses. Their labels encode also reference values for key
physical properties: elasticity and density – see Table I.

Label | Density | Elasticity | Label | Density | Elasticity
---|---|---|---|---|---
V4515 | 45 | 1.5 | V5015 | 50 | 1.5
GV5030 | 50 | 3.0 | GV5040 | 50 | 4.0
N4072 | 40 | 7.2 | NF2140 | 21 | 4.0
T1820 | 18 | 2.0 | T2030 | 20 | 3.0
T3240 | 32 | 4.0 | T2545 | 25 | 4.5
RL3529 | 35 | 2.9 | RL4040 | 40 | 4.0
RL5045 | 50 | 4.5 | RP1725 | 17 | 2.5
RP2440 | 24 | 4.0 | RP27045 | 270 | 4.5
RP30048 | 300 | 4.8 | RP3555 | 35 | 5.5
RP2865 | 28 | 6.5 | RP50080 | 500 | 8.0
TABLE I: Polyurethane foams set. Label – code from manufacturer. Density is in
${kg\cdot m^{-3}}$. Elasticity is the $CV_{40}$ in kPa – “compression stress
value at 40%” [26].

###  III-B Robot grippers and hands

Four different devices were employed, with different morphology and hence
degree of anthropomorphism (industrial type parallel-jaw grippers vs. robot
hands with fingers), different options for controlling the squeezing action
(compression speed and, if available, gripper configuration), and different
type and dimension of feedback. An overview is provided in Fig. 2. Below we
detail technical parameters of every gripper: the action parameters, the
sensory channels, the stopping criterion of the grasping action and the
sampling rate—all these may affect the discrimination of the object being
compressed.

Figure 2: Robot hands and grippers.

An illustration of the data collection process is in the accompanying video.
Fig. 3 shows examples of raw data from individual grippers for action
parameters and objects.

(a) Robotiq 2F-85

(b) QB Soft Hand

(c) OnRobot RG6

(d) Barrett Hand - yellowsponge

(e) Barrett Hand - whitedie

(f) Barrett Hand - yellowcube

Figure 3: Visualisation of raw data for individual grippers. The effect of
action parameters and grasped object. _yellowsponge_ and _yellowcube_ are from
the same material; _yellowcube_ and _whitedie_ have the same size.

####  III-B1 OnRobot RG6

The gripper has a 160 mm stroke and adjustable maximum force threshold. The
gripper fingertips’ surface area is $866\,\mathrm{mm}^{2}$.

Action parameters: Force threshold and closing velocity. We used two force
thresholds: 39 N and 109 N. This choice also affects the gripper closing
speed. The thresholds were chosen such that the closing speed approximately
corresponds to the speeds used on the Barrett Hand (see Section III-B5):
$v1=46$ and $v2=93$ mm/s.

Sensory channels: 2 channels: gripper position (gap between jaws), motor
current.

Stopping criterion: Force threshold reached.

Sampling rate: 15 Hz. An average sample from compression of a single object
contained around 45 points in every channel ($v1$), 41 data points for $v2$.

####  III-B2 Robotiq 2F-85

The gripper has 2 fingers with 85 mm stroke and offers a grip force from 20 N
to 235 N. The fingertips dimensions can be approximated to a rectangle of
$37.5\times 22\,\mathrm{mm}$.

Action parameters: Velocity control was used for the closing. This is quoted
in % of the maximum closing speed. Four different nominal closing speeds were
executed: 0.68% ($v0.0068$, approx. 1.6 mm/s), 20.69% ($v0.21$, 46 mm/s),
63.42% ($v0.63$, 80 mm/s), and 100% ($v1.0$, 131.33 mm/s).

Sensory channels: 2 channels: gripper position (gap between jaws 0–85 cm),
motor current (A).

Stopping criterion: Position does not change anymore (empirical threshold: 8
samples).

Sampling rate: 100 Hz. An average sample from compression of a single object
contained around 760 points in every channel ($v0.0068$), 439 ($v0.21$), 356
($v0.63$), and 319 ($v1.0$) data points, respectively.

####  III-B3 qb SoftHand (Research edition)

This is an anthropomorphic robotic hand with five fingers controlled via a
single electric motor—synchronous opening or closing of all fingers. Thanks to
passive compliance in the joints, it will conform to differently shaped
objects. It can reach a grasp force of 62 N and has a nominal payload of 1.7
kg.

Action parameters: The hand configuration cannot be changed, but one can
choose where to place the object. We considered two actions to minimize
($action1$ or $a1$) and maximize ($action2$ or $a2$) contact with the thumb
(see Fig. 2). The hand is controlled by setting way points to the single
motor. We achieved two closing velocities by setting the closing time to 2.5 s
($v1$) and 1.5 s ($v2$).

Sensory channels: 2 channels: motor position and motor current. Since the
single motor drives all coupled fingers together, these values indirectly (and
ambiguously) code for the position and effort between the fingers and the
object.

Sampling rate: 10 Hz. An average sample from compression of a single object
contained around 28 data points in every channel ($v1$), and 17 data points
for $v2$.

####  III-B4 Barrett Hand

The Barrett Hand (model BH8-282) has three fingers, of which two can rotate
around the base. This design allows numerous configurations of the three
fingers, with the only restriction that the spread joints of Finger 1 and 2
are mirroring each other.

Action parameters: We used two different finger configurations to press the
object in the hand: opposing fingers ($action1$ or $a1$) and lateral
configuration with all fingers on one side pressing against the palm
($action2$ or $a2$) (see Fig. 2). Two objects, blue die and pink die were too
big to be pressed in the lateral finger configuration. Velocity control was
used with two different joint velocities: $v1=0.6$, and $v2=1.2$ rad/s. For
$action2$ we used forward kinematics to calculate an approximation of the
object compression velocity to be 46 and 93 mm/s, respectively.

Sensory channels: The hand is equipped with 96 pressure-sensitive cells (24
per fingertip and 24 on palm) and 3 fingertip torque sensors. In addition, 8
joint angles from encoders were recorded.

Sampling rate: All data was recorded at 200 Hz, but the effective rate of
tactile sensor readings was around 25 Hz. An average sample from compression
of a single object contained around 260 data points in every channel ($v1$),
220 data points for $v2$.

####  III-B5 Compression velocity matching

To the extent that this was possible, the nominal compression velocities (with
empty gripper) were chosen such that they were comparable. The 39 N threshold
on the RG6 results in closing speed of approximately 46 mm/s, which
corresponds to 21% closing speed for the Robotiq 2F-85

##  IV Classification and unsupervised analysis

###  IV-A Time series data processing

As the input for the classifiers we used either: (i) raw sensory data signals;
or (ii) a vector of hand-crafted features. Data collection was organized such
that the number of data points per time series was the same for all trials
(all objects) with a given setting (compression speed). If that was not
possible—trial was terminated when no further compression was possible
(Barrett Hand, Robotiq 2F-85)—the data was padded with zeros from the
beginning.

For this purpose, we made use of features proposed in [9]. These include
minimum, mean, skewness, maximum, kurtosis, median, standard deviation, sum of
values, and amplitude of Hilbert transform. All features are computed for each
sensory channel, concatenated in one multidimensional vector, and normalized
with the zero mean and unit standard deviation across all dimensions of the
feature vector.

###  IV-B Classification

We compare performance of the following classifiers (implementation from
Scikit-learn and PyTorch):

####  IV-B1 $k$-nearest neighbors (kNN)

trained either on raw sensory data or on hand-crafted features. Euclidean
metrics is used and $K$ is selected empirically by evaluating the performance
on the validation set. Each gripper has thus its own optimal $K$ on object set
(objects / foams).

####  IV-B2 Support Vector Machine (SVM)

trained on hand-crafted features. We used grid search to find the optimal
values of hyperparameters (kernel, $C$, and $\gamma$) for a given gripper,
action parameters, and object set.

####  IV-B3 Long-Short Term memory Neural Network (LSTM)

trained on raw data that are padded with zeros to the length of the longest
measurement, keeping the information about the original length. The output of
the last LSTM cell (representing the LSTM-computed features) with nonzero
input is selected. The features are then passed into two linear layers to
compute the output. Object’s category is selected using the softmax layer.
During training, the network weights are modified. Gradient descent (Adam
[27]) is performed to minimize the cross entropy loss. We used grid search to
find the optimal values of the hyperparameters for a given gripper, action
parameters, and object set. The search is through the number of LSTM hidden
layers (2/4), size of the hidden layer (32/64/128/256), learning rate (log.
distribution over [$10^{-5}$, $10^{-3}$]). Batch size was determined
empirically to comply with the available memory resources. This resulted in
336 different LSTM models. The model performing the best on the validation set
was selected (typically networks with 4 LSTM layers and 256 neurons in the
hidden layer).

####  IV-B4 Training, validation, and testing sets

The datasets used were divided into training, validation and testing subsets.
The size of each subset for a given gripper and action parameters is in Table
II (Objects – dataset) or in Table III (Foams – dataset). The dataset sizes
indicate the number of time series—one time series corresponds to a single
object compression. Testing and validation subsets were distinct but had same
size. The number of sensory channels depends on the gripper. The column _all_
is not always the sum of the individual columns—the number of samples was
reduced by using a random selection. The data is available at
https://osf.io/zetg3/?view_only=6f69b0181e394f128a71ce0d811176e6.

qb Soft Hand
---
| a1v1 | a1v2 | a2v1 | a2v2 | all | per object
train | 575 | 604 | 578 | 539 | 1894 | 60+
val/test | 180 | 180 | 180 | 180 | 600 | 20
Barrett Hand
| a1v0.6 | a1v1.2 | a3v0.6 | a3v1.2 | all | per object
train | 666 | 692 | 682 | 692 | 1366 | 74
val/test | 270 | 270 | 270 | 270 | 540 | 30
OnRobot RG6
| 39N | 109N | - | - | all | per object
train | 270 | 270 | - | - | 540 | 30
val/test | 135 | 135 | - | - | 270 | 15
Robotiq 2F-85
| v0.0068 | v0.21 | v0.63 | v1.0 | all | per object
train | 913 | 431 | 373 | 903 | 2620 | 40+
val/test | 135 | 135 | 135 | 135 | 540 | 15
TABLE II: Objects – dataset. Barrett Hand
---
| a1v0.6 | a1v1.2 | a3v0.6 | a3v1.2 | all | per foam
train | 653 | 657 | 658 | 659 | 2627 | 32+
val/test | 300 | 300 | 300 | 300 | 1200 | 15
OnRobot RG6
| 39N | 109N | - | - | all | per foam
train | 600 | 600 | - | - | 1200 | 30
val/test | 300 | 300 | - | - | 600 | 15
Robotiq 2F-85
| v0.0068 | v0.21 | v0.63 | v1.0 | all | per foam
train | 1751 | 628 | 700 | 723 | 3872 | 30+
val/test | 300 | 300 | 300 | 300 | 1200 | 15
TABLE III: Foams – dataset.

###  IV-C Unsupervised analysis of time series data

In order to understand what the dominant sources of variability in the data
are—the different objects explored, the way they are manipulated, the sensory
channels used—unsupervised learning techniques can be employed. For data in
the form of features, standard dimensionality reduction techniques such as
Principal Component Analysis (PCA) can be employed to visualize the data.

##  V Results – classification

###  V-A Ordinary objects set – classification results

Time series from the four different hands/grippers compressing the 9 objects
were used to train and evaluate classifiers (see Section IV-B). An overview
for the different devices is in Fig. 4. The performance of different
classifiers is in rows of every panel. Overall, the choice of classifier does
not play an important role, with k-NN over features performing slightly worse.
Different action parameters—grasp velocity or configuration—are compared in
columns of every panel. Three of the devices—OnRobot RG6, Robotiq 2F-85, and
Barrett Hand—achieve close to perfect performance. The qb Soft Hand is much
less accurate.

(a) OnRobot RG6

(b) Robotiq 2F-85

(c) qb Soft Hand

(d) Barrett Hand

Figure 4: Objects set – classification accuracy. (a) OnRobot RG6. (b) Robotiq
2F-85. (c) QB Soft Hand (d) Barrett Hand.

###  V-B Polyurethane foams set – classification results

In order to pose a bigger challenge, we tested the Robotiq 2F-85, OnRobot RG6,
and the Barrett hand on the set of 20 polyurethane foams (Section III-A2 and
Fig. 1(b)). Note that this is a highly challenging dataset, impossible for
humans to classify correctly. The results are summarized in Fig. 5. There are
now bigger differences between the classifiers, k-NN lagging behind SVM and
LSTM. SVM performs best for the Robotiq 2F-85 gripper (around 70% accuracy)
and the Barrett Hand (95%; 90% with LSTM). OnRobot RG6 achieves best accuracy
(around 75%) with LSTM.

(a) OnRobot RG6

(b) Robotiq 2F-85

(c) Barrett Hand

Figure 5: Polyurethane foams set – classification accuracy. (a) OnRobot RG6
(b) Robotiq 2F-85 (c) Barrett Hand.

###  V-C Effect of action parameters

We studied the effect of speed of object compression. The results are in
columns of Figs. 4 and 5. Overall, there is a trend that faster compression
speeds degrade classification performance. With a fixed sampling rate, slower
speeds generate more data points. However, additional analyses applying up- or
downsampling led to very similar results, suggesting that the mechanical
response of the objects at slower compression speed facilitates
discrimination. Interestingly, for the RG6 gripper, bigger speed was coupled
with higher force threshold and hence greater object deformation. However,
performance was better for a smaller threshold.

On the anthropomorphic hands, different configurations of the hand (Barrett
Hand) or different placement of the object in the hand (qb Soft Hand) were
possible—labeled as actions 1 or 2 (_a1_ or _a2_) in Figs. 4 and 5. The
configuration has an important effect on the qb Soft Hand only, whereby the
configuration _a1_ gives better results.

###  V-D Transfer learning across action parameters

Fig. 6 provides an overview how learning on a dataset with certain action
parameters—compression speed or finger configuration—transfers to data
collected with different parameters. First, if training is performed on all
datasets (last row in individual subplots), the performance on individual
testing sets is as good as that of classifiers trained on that specific
condition (diagonal). Second, for the parallel jaw gripper Robotiq 2F-85 (Fig.
6(a) and 6(b)), there is limited transfer between compression speeds. To some
extent, training on higher speeds partially transfers to lower speeds. On the
Barrett Hand (Fig. 6(e), 6(f), 6(g)), transfer between compression speeds is
much more successful, as the quadrants with the same configuration (_a1_
/_a2_) demonstrate—more so for LSTM compared to SVM.

(a) Robotiq 2F-85 – objects – SVM

(b) Robotiq 2F-85 – foams – SVM

(c) OnRobot RG6 – objects – SVM

(d) OnRobot RG6 – foams – SVM

(e) Barrett Hand – objects – SVM

(f) Barrett Hand – objects – LSTM

(g) Barrett Hand – foams – SVM

Figure 6: Transfer learning – classification accuracy. In individual panes,
rows are marked with training set and columns with testing set. Chance level
performance for the objects set is 11%; for the foams set, it is 5%.

###  V-E Effect of sensory channels (ablation)

We studied which sensory channels are responsible for the classification
performance. Overall, the results are quite intriguing. Even a single sensory
channel is often capable of discriminating the objects. For the parallel jaw
grippers (Fig. 7(a) – Fig. 7(f)), the position channel alone gives very good
performance. The force / current channel is performing worse. For the qb Soft
Hand on the objects set, Fig. 7(c), the pattern is different than for the
parallel jaw grippers—the current sensor is performing better than the
position sensor. On the Barrett Hand, Fig. 7(d) and Fig. 7(g), the tactile
sensors alone provide the best performance, followed closely by position and
effort (fingertip torque) sensors together. Joint position sensors alone
suffice to quite well discriminate the objects due to their different
dimensions Fig. 7(d). On the foams set, Fig. 7(g), position sensors are
importantly more successful in the opposing finger configuration (_a1_). The
fingertip torque sensors alone are overall the worst performing but still
achieving interesting performance.

(a) OnRobot RG6 – objects

(b) Robotiq 2F-85 – objects

(c) qb Soft Hand – objects

(d) Barrett Hand – objects

(e) OnRobot RG6 – foams

(f) Robotiq 2F-85 – foams

(g) Barrett Hand – foams

Figure 7: Effect of sensory channels (ablation) – classification accuracy.
Sensory channels used marked in rows. SVM classifier used in all cases. Chance
level performance for the objects set is 11%; for the foams set, it is 5%.

##  VI Results – Unsupervised analysis

Fig. 8 provides a visualization of the datasets using PCA. Every subplot
visualizes data from one gripper and one object set. If additional action
parameters (velocity, configuration) were available, they are visualized in
the same plot.

(a) OnRobot RG6-objects

(b) OnRobot RG6-foams

(c) Robotiq 2F85-objects

(d) Robotiq 2F85-foams

(e) Barrett Hand-objects

(f) Barrett Hand-foams

(g) qb Soft Hand-objects

(h) 3 grippers-objects

(i) OnRobot RG6 and Robotiq 2F-85 on foams

Figure 8: Visualization using dimensionality reduction (PCA). Colors –
objects. Markers – action parameters.

For the OnRobot RG6 gripper (Fig. 8(a) and Fig. 8(b)), the main clusters are
formed by the action parameters (force threshold / squeezing speed). The
second dimension is due to differences between objects. The softer
objects—yellow cube, yellow sponge, blue cube, and blue die—are clustered
together despite their different dimensions, indicating that stiffness is the
main source of variability.

On the qb Soft Hand, two different configurations and two different velocities
were tried. In Fig. 8(g), the first principal component is dominated by the
velocity parameter, with smaller velocity $v1$ on the left. The second
component is dominated by the objects and in particular their volume.

On the Robotiq 2F85 gripper and the objects set, Fig. 8(c), the first PC is
dominated by the compression velocity, with the two slow speeds grouped
together. The second PC is then taken up by the object characteristics,
dominated by their stiffness. On the foams dataset, Figs. 8(d), the velocity
dominates forming the principal clusters. Note also that the slower speeds
better separate the foams.

On the Barrett Hand and the objects set, the first PC is taken mainly by the
grasping configuration with the opposing fingers ($a1$) on the left side of
Fig. 8(e). Material elasticity seems to dominate PC2 with the softer objects
at the bottom. The compression velocity seems to be the smallest source of
variability. Similar conclusions can be drawn from the visualizations of the
foams dataset (Fig. 8(f)).

The same analysis was performed on the grippers that have two sensory channels
together. In Fig. 8(i), the first principal component separates the two
parallel jaw grippers. Interestingly, the second component separates the
closing speeds, preserving their order. In Fig. 8(h), qb Soft Hand forms a
compact cluster on the center left; OnRobot RG6 on the top. The Robotiq 2F-85
datapoints are most scattered in the bottom part of the plot.

##  VII Discussion and Future Work

We showed that two-dimensional displacement-force data from a single pinch
grasp provides in general sufficient information for discrimination of
deformable objects. Furthermore, this is the case even with uncalibrated
devices and rather low sampling rate (as opposed to e.g. [16]). Ablation
experiments further showed that surprisingly good results can be achieved even
when considering a single sensory channel only. This is in stark contrast with
works that use tactile arrays for object identification [12, 2].

Simple open-loop position or velocity control of squeezing was used and was
found to suffice to discriminate even highly similar objects (polyurethane
foams set). The same approach—but on a professional setup—is used in gauging
material elasticity (e.g., [26]). For situations where one is not handling
delicate objects, such a solution is simpler, quicker, and from a point of
view of material science even better than closed-loop control of the applied
force (e.g., [16]). For single-stream recycling [5], this approach is also the
most efficient, since a single grasp and squeeze suffices for classification
and the object can be immediately sorted to the appropriate container.

Only data from object compression was considered, which is sufficient to gauge
material elasticity. Soft objects like the polyurethane foams, for example,
have an elastic component and a viscous component. Energy absorbed in a
compression/decompression cycle and the dependence of the response on
compression speed reveals also the viscoelastic properties. We experimented
with recording the stress/strain curve also during releasing the grip but no
improvement in classification was achieved. For some devices, technical
limitations did not allow such measurements.

There are several directions of future work. First, not only the compression
but also the decompression phase is relevant for learning about object
properties. Second, the action parameters could be directly fed into the
classifier as additional inputs (discussed in [9]). Finally, a key challenge
and open question is to what extent would transfer learning between the hands
/ grippers be possible. The sensorimotor contingency theory [28] predicts that
long interaction sequences of object exploration make the “object-related
contingencies” stand out from the “apparatus-related contingencies”. That is,
a long enough sequence of sensorimotor interaction with an object may provide
a unique signature of the object’s properties, suppressing the fact that it is
being examined by a hand with very different embodiment.

## References

  * [1] J. Sanchez, J.-A. Corrales, B.-C. Bouzgarrou, and Y. Mezouar, “Robotic manipulation and sensing of deformable objects in domestic and industrial applications: a survey,” _The International Journal of Robotics Research_ , vol. 37, no. 7, pp. 688–716, 2018. [Online]. Available: https://journals.sagepub.com/doi/10.1177/0278364918779698
  * [2] D. Xu, G. E. Loeb, and J. A. Fishel, “Tactile identification of objects using bayesian exploration,” in _2013 IEEE International Conference on Robotics and Automation_. IEEE, 2013, pp. 3056–3061.
  * [3] V. Chu, I. McMahon, L. Riano, C. G. McDonald, Q. He, J. M. Perez-Tejada, M. Arrigo, N. Fitter, J. C. Nappo, T. Darrell _et al._ , “Using robotic exploratory procedures to learn the meaning of haptic adjectives,” in _2013 IEEE International Conference on Robotics and Automation_. IEEE, 2013, pp. 3048–3055.
  * [4] A. J. Spiers, M. V. Liarokapis, B. Calli, and A. M. Dollar, “Single-grasp object classification and feature extraction with simple robot hands and tactile sensors,” _IEEE Transactions on Haptics_ , vol. 9, no. 2, pp. 207–220, 2016.
  * [5] L. Chin, J. Lipton, M. C. Yuen, R. Kramer-Bottiglio, and D. Rus, “Automated recycling separation enabled by soft robotic material classification,” in _2019 2nd IEEE International Conference on Soft Robotics (RoboSoft)_. IEEE, 2019, pp. 102–107.
  * [6] D. Guo, H. Liu, B. Fang, F. Sun, and W. Yang, “Visual affordance guided tactile material recognition for waste recycling,” _IEEE Transactions on Automation Science and Engineering_ , 2021.
  * [7] V. E. Arriola-Rios, P. Guler, F. Ficuciello, D. Kragic, B. Siciliano, and J. L. Wyatt, “Modeling of deformable objects for robotic manipulation: A tutorial and review,” _Fr. Robotics and AI_ , vol. 7, p. 82, 2020.
  * [8] V. Müller and M. Hoffmann, “What is morphological computation? on how the body contributes to cognition and control,” _Artificial Life_ , vol. 23, no. 1, pp. 1–24, 2017.
  * [9] M. Hoffmann, K. Stepanova, and M. Reinstein, “The effect of motor action and different sensory modalities on terrain classification in a quadruped robot running with multiple gaits,” _Robotics and Autonomous Systems_ , vol. 62, no. 12, pp. 1790–1798, 2014.
  * [10] S. Luo, J. Bimbo, R. Dahiya, and H. Liu, “Robotic tactile perception of object properties: A review,” _Mechatronics_ , vol. 48, pp. 54–67, 2017.
  * [11] Q. Li, O. Kroemer, Z. Su, F. F. Veiga, M. Kaboli, and H. J. Ritter, “A review of tactile information: Perception and action through touch,” _IEEE Transactions on Robotics_ , 2020.
  * [12] J. Hoelscher, J. Peters, and T. Hermans, “Evaluation of tactile feature extraction for interactive object recognition,” in _International Conference on Humanoid Robots_. IEEE, 2015, pp. 310–317.
  * [13] M. Kaboli, D. Feng, and G. Cheng, “Active tactile transfer learning for object discrimination in an unstructured environment using multimodal robotic skin,” _International Journal of Humanoid Robotics_ , vol. 15, no. 01, p. 1850001, 2018.
  * [14] M. V. Liarokapis, B. Calli, A. J. Spiers, and A. M. Dollar, “Unplanned, model-free, single grasp object classification with underactuated hands and force sensors,” in _International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2015, pp. 5073–5080. [Online]. Available: https://ieeexplore.ieee.org/abstract/document/7354091
  * [15] L. Scimeca, P. Maiolino, D. Cardin-Catalan, A. P. del Pobil, A. Morales, and F. Iida, “Non-destructive robotic assessment of mango ripeness via multi-point soft haptics,” in _International Conference on Robotics and Automation_. IEEE, 2019, pp. 1821–1826.
  * [16] L. Wang, Q. Li, J. Lam, and Z. Wang, “Tactual recognition of soft objects from deformation cues,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 1, pp. 96–103, 2022.
  * [17] A. Delgado, C. A. Jara, D. Mira, and F. Torres, “A tactile-based grasping strategy for deformable objects’ manipulation and deformability estimation,” in _2015 12th International Conference on Informatics in Control, Automation and Robotics_ , vol. 2. IEEE, 2015, pp. 369–374.
  * [18] K. Hosoda and T. Iwase, “Robust haptic recognition by anthropomorphic bionic hand through dynamic interaction,” in _International Conference on Intelligent Robots and Systems_. IEEE, 2010, pp. 1236–1241. [Online]. Available: https://ieeexplore.ieee.org/abstract/document/5649297
  * [19] M. Kerzel, E. Strahl, C. Gaede, E. Gasanov, and S. Wermter, “Neuro-robotic haptic object classification by active exploration on a novel dataset,” in _Int. Joint Conf. on Neural Networks_. IEEE, 2019, pp. 1–8.
  * [20] M. Madry, L. Bo, D. Kragic, and D. Fox, “St-hmp: Unsupervised spatio-temporal feature learning for tactile data,” in _2014 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2014, pp. 2262–2269.
  * [21] T. Bhattacharjee, J. M. Rehg, and C. C. Kemp, “Haptic classification and recognition of objects using a tactile sensing forearm,” in _International Conference on Intelligent Robots and Systems_. IEEE, 2012, pp. 4090–4097. [Online]. Available: https://ieeexplore.ieee.org/abstract/document/6386142/
  * [22] ——, “Inferring object properties with a tactile-sensing array given varying joint stiffness and velocity,” _International Journal of Humanoid Robotics_ , vol. 15, no. 01, p. 1750024, 2018.
  * [23] S. Decherchi, P. Gastaldo, R. S. Dahiya, M. Valle, and R. Zunino, “Tactile-data classification of contact materials using computational intelligence,” _IEEE Transactions on Robotics_ , vol. 27, no. 3, pp. 635–639, 2011.
  * [24] J. Bednarek, M. Bednarek, P. Kicki, and K. Walas, “Robotic touch: Classification of materials for manipulation and walking,” in _2019 2nd IEEE International Conference on Soft Robotics_. IEEE, 2019, pp. 527–533. [Online]. Available: https://ieeexplore.ieee.org/document/8722819
  * [25] R. L. Klatzky and S. J. Lederman, “Identifying objects from a haptic glance,” _Perception & Psychophysics_, vol. 57, no. 8, pp. 1111–1123, 1995.
  * [26] “Polymeric materials, cellular flexible – Determination of stress-strain characteristics in compression, Part 1: Low-density materials,” 1986. [Online]. Available: https://www.iso.org/obp/ui/#iso:std:iso:3386:-1:ed-2:v1:en
  * [27] D. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _International Conference on Learning Representations_ , 2015.
  * [28] J. K. O’Regan and A. Noe, “A sensorimotor account of vision and visual consciousness,” _Behavioral and Brain Sciences_ , vol. 24, pp. 939–1031, 2001.


# Linking intra- and extra-cellular metabolic domains via neural-network
surrogates for dynamic metabolic control

Sebastián Espinel-Ríos and José L. Avalos Manuscript received October 26,
2023.This work was supported by the U.S. Department of Energy, Office of
Science, Office of Biological and Environmental Research (Award Number DE-
SC0022155).S. Espinel-Ríos and J. L. Avalos are with the Department of
Chemical and Biological Engineering, Princeton University, Princeton, NJ 08544
USA (SER e-mail<EMAIL_ADDRESS>JLA e-mail: javalos@princeton.edu).

###### Abstract

In this study, we aim to optimize biotechnological production by manipulating
intracellular metabolic fluxes in microbial cell factories. Model-based
dynamic optimization is proposed to determine the optimal dynamic trajectories
of the manipulatable intracellular fluxes. A challenge emerges as existing
models are often oversimplified, lacking insights into intracellular
metabolism, or are excessively complex, leading to numerical and
implementation challenges in optimal control (e.g., related to bilevel
optimizations). We propose a solution involving a machine-learning surrogate
derived from steady-state constraint-based metabolic modeling. This surrogate
bridges the gap between manipulatable intracellular fluxes and process
exchange rates. By integrating the surrogate model with simple macro-kinetic
dynamic models, we can develop hybrid machine-learning-supported dynamic
models. Conveniently, the manipulatable intracellular fluxes in these
augmented models can be exploited as dynamic optimization degrees of freedom.
We apply this modeling and optimization strategy to a representative metabolic
network that showcases common challenges in dynamic metabolic control. We also
present an example of cybernetic control to counteract system uncertainties.
Our approach facilitates the in silico evaluation of dynamic metabolic
interventions and can aid in the selection of suitable control and actuation
strategies.

##  1 Introduction

Microbial cell factories are engineered microorganisms that are optimized for
the biotechnological production of valuable metabolites and proteins from
renewable resources. Metabolic engineering involves the rewiring of metabolic
networks to enhance the production of native products or to introduce new
(non-native) metabolic pathways into cells. As a result, these microbial cell
factories can make a wide range of biochemicals, biofuels, biomaterials, and
biopharmaceuticals [1].

Bioprocess efficiency is commonly assessed by the product titer, volumetric
productivity rate, and yield. In conventional metabolic engineering
approaches, the cell’s metabolic steady state is optimized to favor the
product pathway. However, maximizing the product pathway can divert resources
away from biomass synthesis, resulting in a trade-off where increased product
yields lead to decreased biomass yields and reduced volumetric productivities
[2].

Dynamic metabolic control [2] has emerged as a means to address intrinsic
metabolic trade-offs, potentially enhancing the production efficiency of
bioprocesses. It focuses on dynamically adjusting manipulatable metabolic
fluxes to enable transitions between different metabolic states, rather than
maintaining a static metabolic flux distribution. Compared to only exploiting
extracellular/process exchange rates, manipulating intracellular fluxes can
enable more degrees of freedom for advanced metabolic applications, upgrading
the toolbox of bioprocess optimization.

The manipulation of target intracellular metabolic fluxes can be achieved,
e.g., by tuning the expression of enzymes (transcriptional level) or by
modifying enzyme activity (post-transcriptional level) [3]. Control can be
exerted externally, utilizing external signals in a cybernetic approach [4],
or intracellularly, with control mechanisms encoded within the cell through
genetic and molecular circuits [5]. Fig. 1 shows possible control
configurations for dynamic metabolic control.

Figure 1: Representation of dynamic metabolic engineering and possible control
configuration strategies levering, e.g., two manipulatable intracellular
fluxes. System inputs: $\bm{u}\in\mathbb{R}^{n_{u}}$; system measurements:
$\bm{y}\in\mathbb{R}^{n_{y}}$.

A key question arises: what should the optimal dynamic trajectories of the
manipulatable fluxes be to maximize production efficiency? Knowing a priori
the optimal dynamic trajectories of the manipulatable intracellular fluxes can
facilitate, e.g., the determination of whether a cybernetic or intracellular
control approach is more suitable for implementation. Also, it could help to
decide on the most appropriate actuation mechanism, e.g., either at the
transcriptional or post-transcriptional level.

We propose employing model-based dynamic optimization to determine the optimal
dynamic trajectories of the manipulatable intracellular metabolic fluxes that
maximize production efficiency. To do so, we need a suitable dynamic model
that links manipulatable intracellular metabolic fluxes to extracellular
exchange rates111Without loss of generality, here we assume that all products
and substrates are exchanged with the extracellular medium, although one could
also deal with intracellular products. such as substrate uptake, product
excretion, and growth rates.

Unstructured and unsegregated macro-kinetic modeling may be the simplest way
to model dynamic bioprocesses. These models consider biomass as a single
component catalyzing growth and process exchange rates [6]. As such, they do
not explicitly capture information on intracellular metabolism; thus, not
suitable for applications in dynamic metabolic control exploiting
intracellular fluxes.

Constraint-based modeling, such as flux balance analysis (FBA), enables the
prediction of steady-state metabolic flux distributions (cf. [7] for more
details). These (underdetermined) models optimize a cellular objective
function, often assumed to be the maximization of biomass synthesis. The
optimization is constrained by the mass balances, constructed from the
stoichiometric matrix of a metabolic network, which can be derived from gene-
enzyme-reaction associations. Additional constraints related to, e.g.,
thermodynamic feasibility and resource allocation, can also be applied. The
advantage of these models is their simplicity. In their most basic
formulation, they require only the mass balance and irreversibility reaction
constraints. Using constraint-based models, one can estimate maximum
theoretical fluxes or yields under specific conditions and assumed cellular
objectives. However, the steady-state assumption limits the application of
constraint-based modeling in dynamic optimization schemes.

Dynamic versions of constraint-based modeling are available (cf. e.g. [8, 4]).
Integrating them into process optimization, however, turns the task into a
bilevel optimization and requires assumptions on the cell’s dynamic objective
function, e.g., short-term/long-term goals. Furthermore, the numerical
solution of bilevel optimizations often involves game theory-based assumptions
(cf. [9] for more details). One has to decide whether the two optimizations
collaborate (optimistic approach) or conflict (pessimistic approach). The
optimistic approach simplifies the task, allowing the replacement of the inner
optimization with its Karush–Kuhn–Tucker conditions. However, this introduces
non-convexity, even if the original model was convex. Consequently,
constraint-based dynamic models often present substantial challenges in the
context of process optimization.

In this work, we employ neural networks, known for capturing complex
relationships, to bridge the gap between the steady-state solutions of
constraint-based models and macro-kinetic dynamic models. The neural network
effectively creates a direct link between manipulatable intracellular fluxes
and extracellular exchange rates. The resulting model can be regarded as a
hybrid machine-learning supported dynamic model. This methodology simplifies
process optimization by circumventing bilevel optimization schemes while still
enabling the use of manipulatable intracellular fluxes as optimization degrees
of freedom. Our contribution differs from previous approaches [10], where
polynomial-based surrogates were employed to correlate one set of
extracellular fluxes with another, neglecting the exploration and exploitation
of the intracellular domain.

The remainder of this paper is organized as follows. Section 2 introduces our
hybrid machine-learning-supported dynamic metabolic modeling strategy. In
Section 3, we present a model-based dynamic optimization scheme that leverages
manipulatable intracellular fluxes to maximize production efficiency, and we
show how this can be incorporated into cybernetic control approaches. Section
4 serves as a test case, where we consider a representative metabolic network
with two biotechnology-relevant degrees of freedom: one for modulating the
trade-off between growth and product production, and another one for adjusting
the ratio of products of interest.

##  2 Metabolic modeling strategy

We proceed to introduce each of the elements of our proposed metabolic
modeling strategy and show how they are interconnected.

###  2.1 Constraint-based model in steady state

A constraint-based model, often underdetermined, can be expressed as (adapted
from [7]): VF_bio(V), ˙m=SV = 0, V_min ≤V ≤V_max, V_i,min = 0, ∀i ∈I, V_man =
v_man, V_man ⊆V, 0≤c(V).

Here, $\bm{S}\in\mathbb{R}^{n_{m}\times n_{V}}$ is the stoichiometric matrix
associated with the internal metabolites $\bm{m}\in\mathbb{R}^{n_{m}}$ and
metabolic fluxes $\bm{V}\in\mathbb{R}^{n_{V}}$. $\mathbb{I}$ represents the
set of irreversible reactions, hence the lower bound equal to zero. In case of
manipulatable intracellular fluxes, denoted as
$\bm{V_{\mathrm{man}}}\in\mathbb{R}^{n_{\mathrm{man}}}$, they can be
constrained to their corresponding values
$\bm{v_{\mathrm{man}}}\in\mathbb{R}^{n_{\mathrm{man}}}$. Additional (non-
linear) constraints $\bm{c}:\mathbb{R}^{n_{V}}\rightarrow\mathbb{R}^{n_{c}}$
can be integrated to consider factors like resource allocation and
thermodynamics. The assumed cellular objective function is represented by
$F_{\mathrm{bio}}:\mathbb{R}^{n_{V}}\rightarrow\mathbb{R}$, with $\bm{V}$
being the decision variable of the optimization. Note that the model assumes
steady-state conditions of the metabolism (cf. Eq. (2.1)).

###  2.2 Mapping intra- and extra-cellular metabolic fluxes

We categorize metabolic fluxes into intracellular fluxes, denoted as
$\bm{V_{\mathrm{int}}}\in\mathbb{R}^{n_{\mathrm{int}}}$, and extracellular or
exchange fluxes, represented by
$\bm{V_{\mathrm{ext}}}\in\mathbb{R}^{n_{\mathrm{ext}}}$. Thus,
$\bm{V}:=[\bm{V_{\mathrm{int}}}^{\mkern-1.5mu\mathsf{T}},\bm{V_{\mathrm{ext}}}^{\mkern-1.5mu\mathsf{T}}]^{\mkern-1.5mu\mathsf{T}}$.
To systematically explore the impact on the metabolism of varying
intracellular fluxes, we employ a grid search approach. We define sets
$\mathbb{V}_{\mathrm{man},1},\mathbb{V}_{\mathrm{man},2},\ldots,\mathbb{V}_{\mathrm{man},n}$
representing possible values for each of the $n$ manipulatable intracellular
fluxes. The Cartesian product gives us all possible combinations of these flux
values:

$\mathbb{G}=\prod_{i=1}^{n}\mathbb{V}_{\mathrm{man},i}.$ (1)

Each combination is applied as a constraint in Eq. (2.1). Solving the
constraint-based model in (2.1) results in a set of labels
(extracellular/exchange fluxes) and features (manipulatable intracellular
fluxes). A machine-learning approach can be used to learn the mapping between
the intra- and extra-cellular flux domains. Using neural networks, this can be
expressed as:

$\bm{V}_{\mathrm{ext}}=\bm{f_{\mathrm{NN}}}(\bm{V}_{\mathrm{man}},\bm{\Theta}),$
(2)

where
$\bm{f_{\mathrm{NN}}}:\mathbb{R}^{n_{\mathrm{man}}}\times\mathbb{R}^{n_{\Theta}}\rightarrow\mathbb{R}^{n_{\mathrm{ext}}}$
represents the trained neural network and
$\bm{\Theta}\in\mathbb{R}^{n_{\Theta}}$ represents the corresponding
parameters.

Remark. We could optionally also include the remaining (non-manipulatable)
intracellular fluxes
$\bm{V}_{\mathrm{rem}}\in\mathbb{R}^{n_{V}-n_{\mathrm{man}}}$ as labels,
ensuring continued insight into the intracellular metabolic flux distribution:

$\bm{V}_{\mathrm{rem}}=\bm{f_{\mathrm{NN}}}(\bm{V}_{\mathrm{man}},\bm{\Theta}).$
(3)

For simplicity, we consider Eq. (2) in the remainder of this study.

###  2.3 Hybrid dynamic model linked to intracellular fluxes

The neural-network surrogate described above can be integrated into macro-
kinetic dynamic models, effectively creating a hybrid machine-learning-
supported dynamic model. Therein, the process exchange rates are represented
as functions of manipulatable intracellular fluxes. Let
$\bm{z}\in\mathbb{R}^{n_{z}}$ be the extracellular states, including the
biomass dry weight $Bio\in\mathbb{R}$. Without loss of generality, let us
assume a batch process. The dynamics of the system follows:

$\displaystyle\odv{\bm{z}(t)}{t}=Bio(t)\bm{q}(\bm{f_{\mathrm{NN}}}(\bm{V}_{\mathrm{man}}(t),\bm{\Theta}),h(\bm{z}(t),\bm{\theta}_{h})),$
(4a) $\displaystyle
q_{i}:=V_{\mathrm{ext},{i}}(\bm{V}_{\mathrm{man}}(t),\bm{\Theta})h(\bm{z}(t),\bm{\theta}_{h}),$
(4b) $\displaystyle\forall i\in\\{1,2,...,n_{z}\\},$
$\displaystyle\bm{z}(t_{0})=\bm{z_{0}}.$ (4c)

Here,
$\bm{q}:\mathbb{R}^{n_{\mathrm{man}}}\times\mathbb{R}^{n_{\Theta}}\times\mathbb{R}^{n_{z}}\times\mathbb{R}^{n_{\theta_{h}}}\rightarrow\mathbb{R}^{n_{z}}$
describes the biomass-specific exchange reaction rates of the macro-kinetic
model, with $q_{i}\in\bm{q}$. $\bm{\theta}_{h}\in\mathbb{R}^{n_{\theta_{h}}}$
represents the parameters of the function
$h:\mathbb{R}^{n_{z}}\times\mathbb{R}^{n_{\theta_{h}}}\rightarrow\mathbb{R}$,
while $t$ and $t_{0}$ denote the current and initial time, respectively. The
function $h$ accounts for a priori known rate-limiting factors, e.g.,
substrate uptake limitation or product inhibition, often neglected by
constraint-based models222In the absence of (known) rate limitations, one
could optimistically assume $h=1$ when, e.g., the carbon-/energy-source is
available in the medium; otherwise, $h=0$. This simplification can facilitate
the numerical solution of the dynamic system.. As such,
$V_{\mathrm{ext},{i}}$, obtained from Eq. (2), represents maximum theoretical
flux values predicted by the constraint-based model. Henceforth, we will omit
the time-dependency of variables when clear from the context.

##  3 Dynamic metabolic control

We now show how the manipulatable intracellular fluxes can serve as degrees of
freedom for metabolic control via model-based dynamic optimization.

###  3.1 Dynamic optimization

The dynamic optimization problem for maximizing the production efficiency of a
batch process reads: V_man(⋅)—_t_0^t_hJ_p(z,V_man), Eqs. (2),(4a)-(4c), 0
≤g(z,V_man), where
$J_{p}:\mathbb{R}^{n_{z}}\times\mathbb{R}^{n_{\mathrm{man}}}\rightarrow\mathbb{R}$
is the objective function (e.g., volumetric productivity, economic profit,
etc.). The function
$\bm{g}:\mathbb{R}^{n_{z}}\times\mathbb{R}^{n_{\mathrm{man}}}\rightarrow\mathbb{R}^{n_{g}}$
represents possible state and input constraints, addressing economic,
technical, or safety considerations. The degree of freedom of the optimization
problem is the function of manipulatable fluxes $\bm{V_{\text{man}}}(\cdot)$
from $t_{0}$ to the final prediction time $t_{h}$, e.g., the final process
time $t_{f}$ in a batch.

###  3.2 Counteracting system uncertainties: cybernetic approach

Let us consider a cybernetic approach where the manipulatable intracellular
fluxes can be adjusted online using external inputs. Compared to purely in-
cell feedback control approaches, cybernetics offers in general a higher
degree of controllability, adaptability, and tunability. For simplicity, we
assume that the external input $\bm{u_{e}\in\mathbb{R}^{n_{u_{e}}}}$ directly
influences the activity of the metabolic enzymes catalyzing the manipulatable
intracellular fluxes. This matches, e.g., a post-translational control
scenario. Therefore, to achieve a target $\bm{V_{\text{man}}}$, one could
compute the inputs to be applied to the plant via a suitable actuator. This
relationship follows:

$\bm{u_{e}}=\bm{f_{u}}(\bm{V_{\text{man}}},\bm{\theta}_{u}),$ (5)

where
$\bm{f_{u}}:\mathbb{R}^{n_{\mathrm{man}}}\times\mathbb{R}^{n_{\theta_{u}}}\rightarrow\mathbb{R}^{n_{\mathrm{man}}}$.
$\bm{\theta}_{u}\in\mathbb{R}^{n_{\theta_{u}}}$ represents the parameters of
$\bm{f_{u}}$. Here we assume one input per manipulatable intracellular flux.

Remark. In cybergenetic approaches (control at the transcriptional level) [4],
Eq. (5) should be reconsidered to account for possible time-delays arising
from gene transcription and translation processes.

In cybernetic approaches, one could counteract the short-term effects of
model-plant mismatch and disturbances by repeatedly solving the optimization
problem in (3.1). This enables the re-evaluation of the optimal control
problem in real-time, based on the current state of the system, effectively
implementing model predictive control (MPC) [11]. At sampling time $t_{k}$,
for a batch process, this can be formulated as: u(⋅)—_t_k^t_hJ_p(z,V_man),
Eqs. (2),(4a),(4b),(3.1),(5), z(t_k) = z_k.

Note that the decision variable of the optimization changes from the function
$\bm{V_{\text{man}}}(\cdot)$ in (3.1) to the function $\bm{u}(\cdot)$ in (3.2)
since the feedback control is implemented using a cybernetic approach.
Furthermore, $t_{h}$ can move along with the process time, representing a
moving prediction horizon, or remain constant, indicative of a shrinking
prediction horizon. In batch processes with a predetermined final process
time, we often set $t_{h}=t_{f}$, resulting in a shrinking horizon.

##  4 Representative metabolic network

We consider a representative metabolic network (cf. Fig. 2) that captures
common features in metabolic engineering. The network consists of a single
carbon source $A$ and two products of interest $F$ and $G$. We assume that a
specific ratio of the latter products is desirable in the bioprocess for
optimal product formulation, downstream processing, or end-product properties.
For example, this can be the case of biopolymers, where a given ratio of
monomers defines the final polymer characteristics (tensile strength, melting
point, etc.) [12]. In addition, a by-product $D$ is associated with the
biomass formation process due to the necessity of $ATP$ as a biomass
precursor. A core aspect underlined by this network is the trade-off between
biomass formation and the synthesis of products of interest. An illustrative
scenario is when all available carbon flux is channeled through the product
pathway, $V_{4}=V_{\mathrm{ext},{1}}$, leading to the absence of biomass
precursors ($ATP$, $B$, and $C$) and, consequently, no growth. The growth rate
is represented by $V_{\mathrm{ext},{5}}$.

Figure 2: Representative metabolic network with 11 irreversible fluxes
$V_{i}$. Orange arrows indicate extracellular/exchange fluxes. Dotted arrows
denote the pathways leading to the products of interest.
$\bm{m}:=[A,B,C,D,E,F,G,ATP]^{\mkern-1.5mu\mathsf{T}}$,
$\bm{z}:=[A_{\text{ext}},D_{\text{ext}},F_{\text{ext}},G_{\text{ext}},Bio]^{\mkern-1.5mu\mathsf{T}}$.

$V_{4}$ and $V_{6}$ are manipulatable intracellular fluxes,
$\bm{V_{\mathrm{man}}}:=[V_{4},V_{6}]^{\mkern-1.5mu\mathsf{T}}$. Modulating
$V_{4}$ can divert the metabolic flux away from biomass synthesis to favor
product formation. Therefore, $V_{4}$ unlocks a way to modulate the product-
biomass trade-off. On the other hand, the fine-tuning of $V_{6}$ offers a
mechanism to balance the production ratio of $E$ and $G$.

Remark. We express the extracellular states in $\mathrm{mmol/L}$, except for
biomass which is expressed in $\mathrm{g/L}$. The metabolic fluxes are
biomass-specific, expressed in $\mathrm{mmol/g/h}$, except for the growth rate
which is in $\mathrm{g/g/h}$.

###  4.1 Neural-network surrogate of FBA

Following the approach described in Section 2, we define the sets

$\displaystyle\mathbb{V}_{\mathrm{man},1}:=\left\\{0,0.2,0.4,\ldots,V_{\mathrm{ext},{1},\mathrm{max}}\right\\}$
(6a)
$\displaystyle\mathbb{V}_{\mathrm{man},2}:=\left\\{0,0.01V_{4},0.02V_{4},\ldots,0.5V_{4}\right\\},$
(6b)

where $V_{4}\in\mathbb{V}_{\mathrm{man},1}$ and
$V_{6}\in\mathbb{V}_{\mathrm{man},2}$. The maximum element value of
$\mathbb{V}_{\mathrm{man},2}$ is $0.5V_{4}$ to ensure feasibility333Steady-
state mass balances give $V_{\mathrm{ext},{3}}=V_{5}-V_{6}$ and
$V_{5}=V_{4}-V_{6}$. Ensuring $V_{\mathrm{ext},{3}}\geq 0$ (irreversibility)
leads to $V_{5}\geq V_{6}$, hence $V_{6}\leq 0.5V_{4}$.. The Cartesian product
of these sets follows Eq. (1). Each element of $\mathbb{G}$ is a possible pair
$(V_{4},V_{6})$. We solved an FBA problem as described in (2.1), considering
constraints represented by Eqs. (2.1)-(2.1). The assumed cell’s objective
function in (2.1) was set to the maximization of the growth rate, given $0\leq
V_{\mathrm{ext},{1}}\leq V_{\mathrm{ext},{1},\mathrm{max}}$,
$V_{\mathrm{ext},{1},\mathrm{max}}=10$. The explored yield trade-off space
from the FBA solutions is exemplified by Fig. 3.

Figure 3: Yield trade-offs space for production of (a) $F$ against $Bio$ and
(b) $G$ against $Bio$.

We trained a feedforward neural network444One hidden layer with four neurons,
rectified linear unit (ReLU) activation function. 15 % of the data was used
for testing, 80 % of the remaining data was used for training, and 20 % for
validation. to obtain a machine-learning model as in Eq. (2). The coefficient
of determination for the parity plots of the growth rates and product exchange
rates was $R^{2}=1.00$ (not shown), indicating a perfect fit. The latter was
expected as the data generated for the neural-network surrogate was obtained
in silico without any source of system uncertainty.

###  4.2 Hybrid macro-kinetic model

The macro-kinetic dynamic model in batch for the metabolic system described in
Fig. 2 was formulated as in Eqs. (4a)-(4c), with
$\bm{z}:=[A_{\text{ext}},D_{\text{ext}},F_{\text{ext}},G_{\text{ext}},Bio]^{\mkern-1.5mu\mathsf{T}}$.
Note that $q_{i}$ is defined as a function of the manipulatable intracellular
fluxes. Furthermore, we assumed Monod-type substrate limitation, whose
parameters are widely known for several microorganisms and carbon sources. For
the case study, we consider $h$ in Eq. (4b) as:

$h(A_{\text{ext}})=\frac{A_{\text{ext}}}{A_{\text{ext}}+k_{A}},$ (7)

where $k_{A}=0.04\,\mathrm{mmol/L}$ is the assumed substrate affinity
constant.

###  4.3 Dynamic optimization scenarios

We performed several dynamic optimizations using the optimal control problem
defined in (3.1), aiming to maximize the final concentration of the products
of interest, denoted as $J_{p}=F_{\text{ext}}(t_{f})+G_{\text{ext}}(t_{f})$.
The optimization degree of freedom was
$\bm{V_{\text{man}}}(\cdot)|_{t_{0}}^{t_{h}}$ with
$t_{h}=t_{f}=9\,\mathrm{h}$. We furthermore included a constraint demanding a
specific final product ratio:

$r_{G}=\frac{G_{\text{ext}}(t_{f})}{G_{\text{ext}}(t_{f})+F_{\text{ext}}(t_{f})}.$
(8)

The predicted optimal trajectories of manipulatable intracellular fluxes for
six different scenarios, under
$\bm{z_{0}}=[120,0,0,0,0.001]^{\mkern-1.5mu\mathsf{T}}$, are presented in Fig.
4-(a)-(b). These scenarios correspond to S_1: $r_{G}=0.2$, S_2: $r_{G}=0.3$,
S_3: $r_{G}=0.4$, S_4: $r_{G}=0.6$, S_5: $r_{G}=0.7$, and S_6: $r_{G}=0.8$.
The dynamic trajectories of relevant extracellular states for the latter
scenarios are shown in Fig. 4-(c)-(h).

(a) Manipulated $V_{4}$ (b) Manipulated $V_{6}$ (c) S_1 (d) S_2 (e) S_3 (f)
S_4 (g) S_5 (h) S_6

Figure 4: (a)-(b) Predicted optimal trajectories of the manipulatable
intracellular metabolic fluxes for six different scenarios. (c)-(h) Predicted
dynamic profiles of relevant extracellular states for the considered
scenarios.  $A_{\mathrm{ext}}$,  $F_{\mathrm{ext}}$,  $G_{\mathrm{ext}}$,
$\mathrm{Bio}\times 10$.

The optimization predicted for all cases a two-stage fermentation profile,
indicated by the trajectory of $V_{4}$. During about the first 2 h,
$V_{4}\thickapprox 0$ allowed for biomass accumulation in the batch (growth
phase) followed by $V_{4}\thickapprox 10\,\mathrm{mmol/g/h}$, where all flux
was diverted to the product pathway (production phase). Such a two-stage
fermentation addresses the trade-off between biomass growth and product
synthesis toward maximizing product volumetric productivity. $V_{6}$ followed
a more complex trajectory, finely tuned to comply with the demanded final
product ratios.

Remark. We considered piece-wise constant manipulatable intracellular fluxes
of size $0.01\,\mathrm{h}$ to approximate the function
$\bm{V_{\text{man}}}(\cdot)$, otherwise infinite-dimensional. The constraint-
based model was solved in CasADi [13]. The training of the neural network and
the neural-network-supported dynamic optimizations were performed in HILO-MPC
[14].

###  4.4 Metabolic cybernetic case

Now, we examine a scenario where the manipulatable intracellular fluxes are
controlled using a cybernetic approach. Furthermore, we introduce system
uncertainty to the plant by scaling down both $h$ in Eqs. (4a)-(4b) and
$Bio(t_{0})$ by a factor of 0.97. The optimization/controller, on the other
hand, employed the nominal model as described in the previous sections.

Solving a shrinking-horizon MPC with piece-wise constant manipulatable
intracellular fluxes of size $0.5\,\mathrm{h}$, demanding $r_{G}=0.5$, leads
to the results in Fig. 5. The increased sampling time size of
$0.5\,\mathrm{h}$ is more realistic in terms of process monitoring, i.e., for
updating the current state of the plant in a feedback loop. For the sake of
generality, we plot the target manipulated flux values instead of the external
inputs in Fig. 5. We assume that external inputs are available to modulate the
manipulatable fluxes (cf. Eq. (5)).

Overall, MPC finely adjusted the trajectories of the manipulatable
intracellular fluxes to account for the reduced rates and biomass initial
concentration. Both the open-loop optimization and the MPC controller achieved
an $r_{G}=0.50$, as demanded. However, while the open-loop controller reached
$A_{\mathrm{ext}}(t_{f})=45.8\,\mathrm{mmol/L}$, MPC was able to achieve
almost full substrate depletion
$A_{\mathrm{ext}}(t_{f})=0.8\,\mathrm{mmol/L}$, producing about 1.6-fold more
net final product ($F_{\mathrm{ext}}(t_{f})+G_{\mathrm{ext}}(t_{f})$).

Figure 5: (a)-(b) Manipulatable intracellular fluxes and (c)-(f) relevant
state dynamics for the considered cybernetic example in the presence of system
uncertainty. MPC: model predictive control. OLO_mis: open-loop optimization.

##  5 Conclusion

We presented a method that employs machine learning to bridge manipulatable
intracellular metabolic fluxes with process exchange rates. This surrogate
model is trained with data obtained from solutions of a constraint-based
metabolic model. Consequently, a hybrid machine-learning-supported dynamic
model can be formulated, where the exchange process rates depend on the
manipulatable intracellular fluxes. This approach not only enables dynamic
extracellular state predictions given changes in the intracellular metabolism,
but also facilitates dynamic optimization using intracellular fluxes
explicitly as optimization degrees of freedom.

We applied this methodology to a representative metabolic network,
highlighting common challenges in metabolic engineering such as the trade-off
between product and biomass yields, and the balancing of product ratios. Using
this network, we demonstrated our modeling and optimization approach
considering two intracellular metabolic fluxes as degrees of freedom for
dynamic optimization. Additionally, we outlined a cybernetic scheme for real-
time flux adjustments to counteract system uncertainties.

The proposed approach can contribute to the quest for developing advanced
biotechnological processes. It facilitates efficient in silico testing of
diverse dynamic metabolic engineering strategies, actuation mechanisms, and
control approaches. This is especially beneficial when constrained by limited
information or system knowledge.

## References

  * [1] Jae Sung Cho, Gi Bae Kim, Hyunmin Eun, Cheon Woo Moon and Sang Yup Lee  “Designing microbial cell factories for the production of chemicals”  In _JACS Au_ 2.8, 2022, pp. 1781–1799  DOI: 10.1021/jacsau.2c00344
  * [2] Naveen Venayak, Nikolaos Anesiadis, William R Cluett and Radhakrishnan Mahadevan  “Engineering metabolism through dynamic control”  In _Curr Opin Biotechnol_ 34, 2015, pp. 142–152  DOI: 10.1016/j.copbio.2014.12.022
  * [3] Shannon M. Hoffman, Allison Y. Tang and José L. Avalos  “Optogenetics illuminates applications in microbial engineering”  In _Annu Rev Chem Biomol_ 13.1, 2022, pp. 373–403  DOI: 10.1146/annurev-chembioeng-092120-092340
  * [4] Sebastián Espinel-Ríos, Bruno Morabito, Johannes Pohlodek, Katja Bettenbrock, Steffen Klamt and Rolf Findeisen  “Towards a modeling, optimization and predictive control framework for fed-batch metabolic cybergenetics”  In _Biotechnol Bioeng_ , in production
  * [5] Christina V. Dinh and Kristala L.. Prather  “Development of an autonomous and bifunctional quorum-sensing circuit for metabolic flux control in engineered Escherichia coli”  In _Proc Natl Acad Sci USA_ 116.51, 2019, pp. 25562–25568  DOI: 10.1073/pnas.1911144116
  * [6] Sarantos Kyriakopoulos, Kok Siong Ang, Meiyappan Lakshmanan, Zhuangrong Huang, Seongkyu Yoon, Rudiyanto Gunawan and Dong‐Yup Lee  “Kinetic modeling of mammalian cell culture bioprocessing: the quest to advance biomanufacturing”  In _Biotechnol J_ 13.3, 2018, pp. 1700229  DOI: 10.1002/biot.201700229
  * [7] Steffen Klamt, Stefan Müller, Georg Regensburger and Jürgen Zanghellini  “A mathematical framework for yield (vs. rate) optimization in constraint-based modeling and applications in metabolic engineering”  In _Metab Eng_ 47, 2018, pp. 153–169  DOI: 10.1016/j.ymben.2018.02.001
  * [8] Steffen Waldherr, Diego A. Oyarzún and Alexander Bockmayr  “Dynamic optimization of metabolic networks coupled with gene expression”  In _J Theor Biol_ 365, 2015, pp. 469–485
  * [9] Stephan Dempe  “Bilevel optimization: theory, algorithms, applications and a bibliography”  In _Bilevel Optimization_ 161 Cham: Springer International Publishing, 2020, pp. 581–672
  * [10] Rafael D. De Oliveira, Matheus N. Guedes, José Matias and Galo A.. Le Roux  “Nonlinear predictive control of a bioreactor by surrogate model approximation of flux balance analysis”  In _Ind Eng Chem Res_ 60.40, 2021, pp. 14464–14475  DOI: 10.1021/acs.iecr.1c01242
  * [11] James Rawlings, David Mayne and Moritz Diehl  “Model predictive control: theory, computation and design”  Santa Barbara: Nob Hill Publishing, LLC, 2020
  * [12] T. Volova, E. Kiselev, I. Nemtsev, A. Lukyanenko, A. Sukovatyi, A. Kuzmin, G. Ryltseva and E. Shishatskaya  “Properties of degradable polyhydroxyalkanoates with different monomer compositions”  In _Int J Biol Macromol_ 182, 2021, pp. 98–114
  * [13] Joel A.. Andersson, Joris Gillis, Greg Horn, James B. Rawlings and Moritz Diehl  “CasADi: a software framework for nonlinear optimization and optimal control”  In _Math Program Comput_ 11.1, 2019, pp. 1–36
  * [14] Johannes Pohlodek, Bruno Morabito, Christian Schlauch, Pablo Zometa and Rolf Findeisen  “Flexible development and evaluation of machine-learning-supported optimal control and estimation methods via HILO-MPC”  In _arXiv:2203.13671_ , 2022