# TransMA: an explainable multi-modal deep learning model for predicting
properties of ionizable lipid nanoparticles in mRNA delivery

Kun Wu11footnotemark: 1<EMAIL_ADDRESS>Zixu Wang22footnotemark: 2
<EMAIL_ADDRESS>Xiulong Yang<EMAIL_ADDRESS>Yangyang Chen
<EMAIL_ADDRESS>Zhenqi Han<EMAIL_ADDRESS>Jialu Zhang
<EMAIL_ADDRESS>Lizhuang Liu<EMAIL_ADDRESS>
###### Abstract

As the primary mRNA delivery vehicles, ionizable lipid nanoparticles (LNPs)
exhibit excellent safety, high transfection efficiency, and strong immune
response induction. However, the screening process for LNPs is time-consuming
and costly. To expedite the identification of high-transfection-efficiency
mRNA drug delivery systems, we propose an explainable LNPs transfection
efficiency prediction model, called TransMA. TransMA employs a multi-modal
molecular structure fusion architecture, wherein the fine-grained atomic
spatial relationship extractor named molecule 3D Transformer captures three-
dimensional spatial features of the molecule, and the coarse-grained atomic
sequence extractor named molecule Mamba captures one-dimensional molecular
features. We design the mol-attention mechanism block, enabling it to align
coarse and fine-grained atomic features and captures relationships between
atomic spatial and sequential structures. TransMA achieves state-of-the-art
performance in predicting transfection efficiency using the scaffold and cliff
data splitting methods on the current largest LNPs dataset, including Hela and
RAW cell lines. Moreover, we find that TransMA captures the relationship
between subtle structural changes and significant transfection efficiency
variations, providing valuable insights for LNPs design. Additionally,
TransMA’s predictions on external transfection efficiency data maintain a
consistent order with actual transfection efficiencies, demonstrating its
robust generalization capability. The code, model and data are made publicly
available at https://github.com/wklix/TransMA/tree/master. We hope that high-
accuracy transfection prediction models in the future can aid in LNPs design
and initial screening, thereby assisting in accelerating the mRNA design
process.

###### keywords:

ionizable lipid nanoparticles , multi-modal molecular structure ,
interpretability , transfection cliffs

††journal: Nuclear Physics B

[label1]organization=Shanghai Advanced Research Institute, Chinese Academy of
Sciences, city=Shanghai, postcode=201210, country=China
[label2]organization=Department of Computer Science, University of Tsukuba,
city=Tsukuba, postcode=3058577, country=Japan [label3]organization=Hubei
Provincial Key Laboratory of Artificial Intelligence and Smart Learning,
Central China Normal University, city=Wuhan, postcode=430079, country=China
[label4]organization=University of Chinese Academy of Sciences, city=Beijing,
postcode=100049, country=China [label5]organization=School of Computer
Science, Central China Normal University, city=Wuhan, postcode=430079,
country=China co-firstco-firstfootnotetext: These authors contributed equally
to this work.

An explainable LNPs transfection efficiency prediction model called TransMA.

TransMA employs a multi-modal molecular structure fusion architecture.

TransMA reveals the atomic-level structure-transfection relationships.

TransMA achieves state-of-the-art performance on the largest current LNPs
dataset.

##  1 Introduction

mRNA-based technologies hold promise for therapeutic applications in fields
such as viral vaccines, protein replacement therapies, cancer immunotherapies,
genome editing[1, 2, 3, 4, 5, 6], and have the potential to reshape the
landscape of life science research and medicine[7, 8, 9]. However, achieving
targeted delivery and intracellular release from endosomes remains challenging
for mRNA delivery systems, emphasizing the critical demand for safe and
effective mRNA delivery materials. Specifically, lipid nanoparticles (LNPs)
have undergone extensive investigation and have successfully transitioned into
clinical use for the delivery of mRNA[10, 11, 12]. For example, mRNA1273 and
BNT162b21 utilize lipid nanoparticles for the delivery of antigen mRNA[13,
14]. LNPs consist of four components: ionizable lipid, phospholipid,
cholesterol, and PEGylated lipid. Among these, the ionizable lipid component
holds the highest molar ratio, which determines the formulation’s delivery
efficiency and stability, serving as the core structure of LNPs[15].
Therefore, the key to selecting for efficient LNPs lies in the selection of
the ionizable lipid[16, 17].

Although component chemistry methods based on three-component reactions
(3-CR)[18, 19] and four-component reactions (4-CR)[20] can synthesize a number
of ionizable lipids, manually testing the transfection efficiency of each
synthesized lipid is time-consuming and costly[21, 22]. Several works have
shown that artificial intelligence represented by machine learning and deep
learning can achieve automatic prediction of LNPs transfection efficiency. For
instance, Ding et al.[9] employed four machine learning methods—support vector
machine, random forest, eXtreme gradient boosting, and multilayer
perceptron—to classify the transfection efficiency of 572 lipid nanoparticles
(LNPs) into two categories, achieving a classification accuracy of 98$\%$.
AGILE [23] is a prediction platform for LNPs transfection efficiency based on
graph convolutional neural networks and pre-training, achieving a mean squared
error (MSE) of approximately 6 in predicting LNPs transfection efficiency.
TransLNP[24] which is a model based on Transformer architecture and data
balancing achieved a MSE of approximately 5 on the AGILE dataset.

Despite significant progress in previous works, it must be acknowledged that
there are notable limitations in achieving automatic prediction of LNPs
transfection efficiency based on machine learning and deep learning. These
limitations hinder the accuracy and generalizability of the models. These
limitations include:(1)The lack of multi-modal information interaction.
Previous works only extract single-modal information from molecules, which
limits their ability to capture Quantitative Structure-Activity Relationships
(QSAR) between ionizable lipids and transfection efficiency. For instance,
chemberta [25] relies on one-dimensional representations to predict molecular
properties. Integrating multi-modal information from molecules is crucial for
further improving prediction accuracy. (2)Previous works lacks
interpretability. Although previous works can predict transfection efficiency
accurately, it is unclear which atoms play a key role in this process. This
lack of interpretability restricts deep understanding of the prediction
process and hampers efforts to further optimize the model for improved
accuracy. Therefore, enhancing model interpretability while maintaining
prediction accuracy poses a significant challenge in current research.
(3)Limited attention to the impact of transfection cliffs. Pairs of molecules
that are structurally very similar but have significantly different
transfection efficiencies—known as transfection cliffs—capture the knowledge
hidden in the elusive structure-property relationships but have received
limited attention.

To address the challenges, we introduce an explainable LNPs transfection
efficiency prediction model called TransMA, which adopts a multi-modal
molecular structure fusion architecture. TransMA integrates 3D geometric
information and 1D atomic sequence information of molecules to predict LNPs
transfection efficiency. A self-attention mechanism SE(3) Transformer
architecture named molecule 3D Transformer is employed to extract 3D geometric
information. Molecule 3D Transformer is pre-trained by reconstructing atomic
3D coordinates and masked atom prediction. To extract atomic sequence
information, the state space model designated as molecule Mamba is presented
to obtain coarse-grained atomic sequence information. By constructing multi-
level molecular structure data pairs, the model jointly learns the multi-
dimensional structural features of molecules. Specifically, we design the mol-
attention mechanism block to align and concatenate the fusion features. The
fusion features include atomic sequence distribution information, atomic
coordinates information, relative atomic positions information, and types of
bonds between atoms information. Compared to state-of-the-art molecular graph
convolutional networks and Transformer models, TransMA demonstrates the best
performance, reducing MSE by 43$\%$ compared to the baseline model AGILE and
by 35$\%$ compared to TransLNP. Additionally, TransMA possesses
interpretability to better understand the mapping relationship between
molecular structure and biochemical properties[26, 27]. The mol-attention
mechanism block not only integrates multi-modal molecular features but also
reveal critical sites that influence transfection efficiency. Transfection
cliffs represent pairs of molecules that are structurally similar but exhibit
substantial differences in transfection efficiency. Extracting these
transfection cliff data can help us understand the elusive transfection
efficiency. A total of 4267 and 2104 transfection cliffs are screened from
Hela and RAW 264.7 cell lines, respectively, constituting 68$\%$ and 81$\%$ of
all data, indicating the widespread occurrence of transfection cliffs. Our
findings is that atoms with high attention scores calculated by the mol-
attention mechanism block corresponded to key atoms in transfection cliff
pairs. This finding indicates that significant differences in transfection
efficiency caused by subtle structural changes. To test the generalization
capability of TransMA, the results demonstrate that TransMA’s predicted values
maintain a consistent ranking with the actual transfection efficiency values
without training on the external dataset.

In this study, our primary contributions can be summarized as follows:

  1. 1.

We propose an explainable LNPs transfection efficiency prediction model, named
TransMA, which adopts a multi-modal molecular structure fusion architecture.
The fine-grained atomic spatial relationship extractor named molecule 3D
Transformer captures molecular 3D geometric features, while the coarse-grained
atomic sequence extractor named molecule Mamba captures 1D SMILES
representation molecular features.

  2. 2.

TransMA achieves state-of-the-art performance on the current largest LNPs
dataset, showing significant improvement across various metrics compared to
previous methods.

  3. 3.

TransMA demonstrates strong interpretability. We construct 4267 and 2104 pairs
of transfection cliffs in Hela and RAW 264.7 cell lines to validate the
model’s identification of key molecular structures. While validating TransMA’s
prediction results, we find that it can detect substantial transfection
efficiency differences caused by minor structural variations (Figure 7).

  4. 4.

We collect independent external test data from several sources. TransMA’s
predicted values consistently maintain the same order as the actual
transfection efficiencies, highlighting its robust generalization capacity.

##  2 Materials and methods

###  2.1 The method of TransMA to predict LNPs transfection efficiency

Deep learning is a powerful tool for accelerating research in the molecular
field[28, 29]. We develop an explainable deep learning approach to predict the
transfection efficiency of LNPs. This model not only achieves high-precision
predictions but also elucidates the key molecular structures that the model
focuses on. The overall architecture of TransMA is illustrated in figure 1.
The TransMA framework consists of three components: molecule 3D Transformer,
molecule Mamba, and mol-attention mechanism block.

Figure 1: The model takes as input multimodal structural information of
ionizable lipids, including three-dimensional structural details: atomic type
sequences, three-dimensional coordinates, distance matrices, bond type
matrices, and one-dimensional SMILES representation molecular structural
information. These two types of structural information undergo feature
extraction separately using self-attention mechanism models and spatial state
models. The resulting features are fused and visualized through a mol-
attention mechanism to reveal the atoms that have a significant impact on
transfection efficiency prediction.

####  2.1.1 Molecule 3D Transformer

The molecule 3D Transformer is employed to extract the three-dimensional
structural information features of molecules. These three-dimensional
structural details include the atomic type sequences, atomic three-dimensional
coordinates, Euclidean distance matrices between atoms, and matrices
representing the types of bonds between atoms in ionizable lipid molecules.
TransLNP adopts a pre-training followed by fine-tuning training approach shown
in figure 2A. Despite the scarcity of labeled data with LNPs properties in the
LNPs field, we can leverage various molecular structure information from other
small molecules to train the model. The pre-training task utilizes masked
language modeling: randomly masking 15$\%$ of the atoms and introducing mask
tokens. Subsequently, an embedding layer is used to map the atomic sequences
to representations of atomic types, providing each atom with semantic
information. Equation 1 illustrates the process of self-attention mechanism
addressing interactions between atoms:

$attention=\text{softmax}(\frac{QK}{\sqrt{d}}+bias)V,$ (1)

The masked atomic types are linearly represented as Queries (Q), Keys (K), and
Values (V) in the self-attention module. $\sqrt{d}$ represents the dimensions
of vectors Q, K, and V. The attention bias in the self-attention mechanism is
based on representations of atomic distances and bond types.

After considering the relationships between atoms, uniform noise is introduced
and added to the 3D coordinates of the molecule. This process of coordinate
updating effectively incorporates the noise into the molecular representation.
In the pre-training phase, atomic type heads, atomic coordinate heads, and
atomic pairwise Euclidean distance heads are employed to predict masked atomic
types, coordinates, and relative distances, respectively. Smooth L1 is used as
the loss function for predicting Euclidean distances and masked atomic
coordinates, while cross-entropy is used for predicting atomic types. During
fine-tuning, the pre-trained model is loaded, following the data processing
procedures from pre-training. Finally, molecular prediction heads are used for
transfection efficiency prediction.

####  2.1.2 Molecule Mamba

Mamba is a linear-time sequence modeling method based on selective state
spaces[30]. Variants of Mamba have demonstrated outstanding performance in
natural language processing, image processing, remote sensing, and audio
domains[31, 32, 33, 34]. State Space Sequence Models (SSMs)[35] form the
theoretical foundation of Mamba, capable of mapping a one-dimensional function
or sequence $u(t)$ to $y(t)\in\mathbb{R}$ through a hidden state $x(t)$. This
mapping can be represented by the following linear ordinary differential
equation shown in equation 2:

$\dot{x}(t)=Ax(t)+Bu(t),\quad y(t)=Cx(t),$ (2)

where the state matrix $A\in\mathbb{R^{N\times N}}$ serves as the evolution
parameter and $B\in\mathbb{R^{N\times 1}}$, $C\in\mathbb{R^{1\times N}}$ act
as projection parameters, representing the implicit latent state.

Due to its combination of recurrent neural networks (RNNs) and convolutional
neural networks (CNNs), Mamba has advantages in handling long sequence data.
The SMILES representation of ionizable lipids can be viewed as a long sequence
for processing. Therefore, we first apply Mamba to the field of molecular
property prediction, using one-dimensional SMILES representation of ionizable
lipids as sequence inputs. Figure 2B illustrates the overall process of
Molecule Mamba. Molecule Mamba involves four parameters
$(Delta(\Delta),A,B,C)$. Firstly, the feature dimension of the SMILES
representation molecule sequence is doubled through linear projection. After
the split operation, the SMILES representation molecule is represented in two
parts. One part undergoes one-dimensional convolution and Sliu activation
function to extract molecular features. After linear projection and split
operation, parameters (Delta, B, C) are obtained. Delta is updated through
softplus. At the same time, parameters (A, D) are initialized, where A and D
are independent of the input. $A$ is discretized using zero-order hold (ZOH)
discretization, while $B$ is discretized using a simplified Euler
discretization shown in equation 3.

$\overline{A}=\exp(\Delta A),\overline{B}=x+\Delta\cdot f(x,u),$ (3)

Combining equation 2 with the scanning operation, the sequential processing of
elements in the SMILES representation molecule involves updating atoms at each
time step based on the cumulative effect of the previous atom and the current
input, effectively propagating information across the entire sequence.
Equation 4 describes this process. Finally, the output is the multiplication
of another portion of the input with the output of the scanning operation.

$\dot{x}(t)=\overline{A}x(t)+\overline{B}u(t),\quad y(t)=Cx(t),$ (4)

Figure 2: Left:The architecture of the molecule 3D Transformer adopts a pre-
training and fine-tuning approach. The pre-training task involves masked
language modeling and noise coordinate prediction. Fine-tuning is conducted
using molecular prediction heads for transfection efficiency prediction.
Right: The input of Molecule Mamba is the SMILES representation of ionizable
lipids. Through the selective scan operation and the discretized SSM state
equation, high-dimensional features of the molecule are extracted.

####  2.1.3 Mol-attention mechanism block

The purpose of the mol-attention mechanism block is to integrate the output
features of the molecule 3D Transformer and molecule Mamba, and introduce an
attention mechanism to explain the model’s attention to each atom. The output
feature $z_{1}$ of the molecule 3D Transformer has a size of $[n,512]$, where
$n$ represents the number of atoms in the ionizable lipid. Meanwhile, the
output feature $z_{2}$ of molecule Mamba has a size of $[m,512]$, where $m$
denotes the sequence length of the SMILES representation molecule. Apart from
containing atomic features, $z_{2}$ also includes features of characters in
the SMILES representation molecule. Before concatenating the features, all
atomic features $z_{2}^{\prime}$ need to be extracted from $z_{2}$ to match
the size of $z_{1}$. Therefore each atomic feature in $z_{2}^{\prime}$
corresponds to one atomic feature in $z_{1}$, ensuring correct alignment
during feature fusion. Like attention mechanisms in the image domain[36, 37,
38], the mol-attention mechanism block aims to weight each atom, compressing
the fused features containing global information of size $n\times 1024$
directly into a $n\times 1$ feature vector $z$. The first fully connected
layer $W_{1}$ compresses the 1024 concatenated features into $1024/r$ channels
to reduce computational complexity. Following a ReLU non-linear activation
layer $\delta$, the second fully connected layer $W_{2}$ produces weights
$attention\ scores$ through a Sigmoid activation $\sigma$. The dimension of
$attention\ scores$ obtained is $n\times 1$, which is shown as equation 5:

$z=\text{concat}(z_{1},z_{2}^{\prime}),attention\
scores=\sigma(g(z,W))=\sigma(W_{2}\delta(W_{1}z)),$ (5)

where $r$ denotes the compression ratio. Finally, the scale operation and a
regression layer are applied to predict the transfection efficiency, which is
shown as equation 6:

$transfection\ efficiency=\tanh(W_{2}(\tanh(W_{1}(z\cdot attention\
scores)))),$ (6)

where $\tanh$ denotes the tanh activation function.

###  2.2 Loss function

The loss function of the TransMA model combines the mean squared error (MSE)
loss function and the triplet loss function, constructing a hybrid loss
function. Equation 7 is shown as:

$hybrid\\_loss=\frac{1}{N}\sum_{i=1}^{N}(y_{i}-\hat{y}_{i})^{2}+\beta\times\mathcal{L}_{\text{triplet}}(z_{1},z_{2}^{\prime}),$
(7)

where $\beta$ is a weight parameter, $N$ is the number of samples, $y_{i}$ and
$\hat{y}_{i}$ respectively represent the true values and the predicted values,
and the definitions of $z_{1}$ and $z_{2}^{\prime}$ are consistent with
equation 5. Equation 8 describes the triplet loss:

$\begin{split}\mathcal{L}_{\text{triplet}}=\frac{1}{\text{num\\_positive\\_losses}+\epsilon}\sum_{i=1}^{2batch}\sum_{j=1}^{2batch}\sum_{k=1}^{2batch}\\\
\max(d(z_{1i},z_{1j})-d(z_{1i},z_{2k})+\text{margin},0)\cdot\text{mask}_{ijk},\end{split}$
(8)

where $batch$ represents the batch size, $z_{1i}$ represents the feature
extracted from the molecule 3D Transformer for the $i$th sample,
$z_{2k}^{\prime}$ denotes the feature extracted from the molecule Mamba for
the $k$th sample, $d(\cdot)$ denotes the Euclidean distance,
num_positive_losses is the number of positive losses, $\epsilon$ is a small
number used for numerical stability, $margin$ is a hyperparameter controlling
the margin size for effective triplets. and $\text{mask}_{ijk}$ is a triplet
mask used to filter out invalid triplets.

The Triplet Loss first computes the Euclidean distance between $z_{1}$ and
$z_{2}$. For each sample, it constructs all possible triplets and computes the
loss value for each triplet. This loss value is calculated as the difference
between the distance from the anchor sample to the positive sample and the
distance from the anchor sample to the negative sample, with a margin added,
and then taking the maximum value. Subsequently, invalid triplets are filtered
out using a triplet mask, and the average loss value is computed as the final
result of the triplet Loss. Triplet loss maps the molecular 3D structural
features extracted by the molecule 3D Transformer and molecule Mamba, along
with the one-dimensional sequence features, to a unified embedding space. By
learning the similarity of molecules in this embedding space, the model can
better generalize to unseen molecules.

##  3 Experimental Section

###  3.1 Dataset

The dataset for TransMA predicting LNPs transfection efficiency comprises two
main components. One part is utilized for pretraining the molecule 3D
Transformer, while the other part is used for training the molecule Mamba and
fine-tuning the molecule 3D Transformer. The dataset used for pretraining the
molecule 3D Transformer is sourced from Unimol[43], which comprises 19 million
diverse molecular structures333https://bioos-hermite-beijing.tos-cn-
beijing.volces.com/unimol_data/pretrain/ligands.tar.gz. The dataset utilized
for training the molecule Mamba and fine-tuning the molecule 3D Transformer is
obtained from AGILE444https://github.com/bowang-lab/AGILE, consisting of 1200
ionizable lipid SMILES representation molecules along with their corresponding
transfection efficiency data in Hela and RAW 264.7 cell lines.

Ionizable lipids comprises 20 head groups, 12 carbon chains with ester bonds,
and 5 carbon chains with isocyanide head groups. High-throughput synthesis of
LNPs entails mixing an aqueous phase containing mRNA and an ethanol phase
containing lipids using the OT-2 pipetting robot. The aqueous phase, prepared
in pH 4.0 10 mM sodium citrate buffer, contains firefly luciferase mRNA, Cre
recombinase mRNA, or EGFP-mRNA. Meanwhile, the ethanol phase comprises a
mixture of 1200 ionizable lipids and fixed proportions of helper phospholipids
(DOTAP, DOPE, cholesterol, and C14-PEG 2000), with a lipid-to-mRNA weight
ratio of 10:1. Transfection efficiency data were obtained by applying these
1200 mRNA-LNPs to Hela and RAW 264.7 cell lines for in vitro transfection
experiments.

###  3.2 Experimental processing

The tasks of molecule 3D Transformer pretraining include using random
positions as corrupted input 3D positions and training the model to predict
the correct positions, employing different heads to predict distances between
corrupted atom pairs, the correct coordinates of corrupted atoms, and to mask
and predict the atom types of corrupted atoms. Atomic type prediction employs
the cross-entropy loss function with a weight of 1. Prediction related to
atomic coordinates and interatomic distances utilizes the smooth L1 loss
function with weights set to 5 and 10 respectively. Pretraining parameters are
set as follows: model layers are 15, batch size is 128, atom types are 30,
learning rate is adjusted using linear decay, and the optimizer is Adam with
$\epsilon$. Parameters for molecule 3D Transformer fine-tuning are as follows:
model layers are 16, batch size is 4, epochs are set to 200, and the initial
learning rate is set to 1e-5.

Molecule Mamba first tokenizes 1200 ionizable lipid SMILES molecular
representation. The tokenizer used is ChemBERTa-77M-MTR from Deepchem. The
purpose of tokenization is to segment SMILES molecule repren into tokens and
convert these tokens into numerical representations that the model can
process. The parameters of the Molecule Mamba model are set as follows: the
feature dimension is 512, the number of layers is 2, the vocabulary size is
100, and the training process parameters are the same as molecule 3D
Transformer fine-tuning. In TransMA’s hybrid loss function, the parameter
$\beta$ is set to 6 during the scaffpld data splitting and 3 during the cliff
data splitting.

###  3.3 Comparison with representative deep learning-based molecular property
prediction models

To demonstrate the superiority of the proposed method, compare the predictive
transfection efficiency accuracy of TransMA with advanced molecular property
prediction models based on graph convolutional neural networks and Transformer
on the LNPs dataset. The compared five models include AGILE, large-scale self-
supervised pretraining for molecular property prediction model named
ChemBERTa[25], self-supervised graph neural network framework named
MolCLR[39], and geometry enhanced molecular representation learning model
named GEM[40], TransLNP.

The LNPs dataset includes 1200 ionizable lipid SMILES molecules and
transfection efficiency data obtained from Hela and RAW 264.7 cell lines. In
model comparison, the LNPs dataset is divided into cliff[41] and scaffold data
splitting approaches. Figure 3 represents LNPs dataset distributions under two
data splitting methods in Hela and RAW 264.7 cell lines. Scaffold splitting
divides LNPs dataset into training set, validation set and testing set as
shown in figure 3A and figure 3D. Cliff splitting utilizes molecular ECFP
descriptors to represent their structures, employing algorithms like spectral
clustering to partition molecules into five clusters. Cliff splitting divides
LNPs dataset into training set, validation set and testing set as shown in
figure 3B and figure 3E. For each cluster, a stratified sampling strategy is
employed to allocate 10$\%$ of the molecules to the testing set and other
molecules to the training set and validation set in an 8:2 ratio. Cliff
splitting ensures that training and test sets contain the proportion of cliff
molecule pairs in the training and testing sets while preserving information
about molecular structural similarity. Scaffold splitting involves grouping
molecules based on their core structural scaffolds, ensuring that both
training and test sets contain molecules from diverse scaffold structures[42].
Cliff splitting considers the similarity of data, while scaffold splitting
considers the diversity of data. Therefore, combining both cliff and scaffold
data splitting methods can comprehensively validate the model’s generalization
capability.

Figure 3:  The distribution of transfection efficiency dataset under two data
splitting methods in Hela and RAW 264.7 cell lines is as follows: the outer
ring of the pie chart represents the training set distribution, the middle
ring represents the validation set distribution, and the inner ring represents
the test set distribution.(A) The distribution of transfection efficiency
dataset under the scaffold-based train-val-test split in Hela cell line . (B)
The distribution of transfection efficiency dataset under the cliff-based
train-val-test split in Hela cell line. (C) The distribution of transfection
efficiency dataset in Hela cell line. (D) The distribution of transfection
efficiency dataset under the scaffold-based train-val-test split in RAW 264.7
cell line. (E) The distribution of transfection efficiency dataset under the
cliff-based train-val-test split in RAW 264.7 cell line. (F) The distribution
of transfection efficiency dataset in RAW 264.7 cell line.

Table 1 presents the results of TransMA and the compared models in predicting
transfection efficiency on the LNPs dataset. Evaluation metrics for assessing
prediction performance include mean squared error (MSE), mean absolute error
(MAE), the coefficient of determination R2, and the Pearson correlation
coefficient (PCC). The results demonstrate that TransMA exhibits superior
performance compared to the five other models. Figure 4 shows the boxplots of
mean squared error (MSE) for six models in Hela and RAW 264.7 cell lines. It
can be observed that TransMA exhibits the lowest MSE in predicting
transfection efficiency on both cell lines, indicating the smallest
fluctuation in prediction errors. Figure 5 show TransMA’s capability to
extract multi-modal fused molecular features under the scaffold and cliff data
splitting methods in Hela and RAW 264.7 cell lines. It can be observed that
the molecular features extracted by molecule 3D Transformer and molecule Mamba
exhibit discrete states in different transfection efficiency intervals. After
feature fusion by TransMA, the molecular features of TransMA show
localization, meaning that the distribution of features for the same
transfection efficiency interval is close. Figure 5 indicates that TransMA has
successfully established a mapping relationship between the ionizable lipid
structure and transfection efficiency.

Figure 4:  Comparison of TransMA performance in predicting transfection efficiency with five different models under the scaffold and cliff data splitting methods for Hela and RAW 264.7 cell lines. (A) The box plots comparing TransMA with five models in Hela cell line. (B) The box plots comparing TransMA with five models in RAW 264.7 cell line. Table 1:  The results of predicting LNPs transfection efficiency in both Hela and Raw 264.7 cell lines by six models under scaffold and cliff data splitting methods are presented. Prediction accuracy is assessed using MSE (smaller values indicate better performance), MAE(smaller values indicate better performance), R2 (larger values indicate better performance), and PCC (Pearson correlation coefficient, larger values indicate better performance). ”BERTa” refers to the chemBERTa model, while ”RAW” refers to Raw 264.7 cells. The best results are highlighted in bold.  Splitting |  | TransMA | AGILE | BERTa | MolCLR | GEM | TransLNP
---|---|---|---|---|---|---|---
| Hela | RAW | Hela | RAW | Hela | RAW | Hela | RAW | Hela | RAW | Hela | RAW
Scaffold | MSE$\downarrow$ | 3.64 | 1.63 | 6.38 | 2.07 | 4.38 | 2.11 | 6.18 | 2.00 | 5.89 | 2.24 | 5.58 | 1.88
MAE$\downarrow$ | 1.57 | 0.98 | 1.92 | 1.11 | 1.70 | 1.13 | 1.98 | 1.09 | 1.92 | 1.16 | 1.92 | 1.06
$R^{2}$ $\uparrow$ | 0.49 | 0.23 | 0.11 | 0.02 | 0.39 | 0.01 | 0.12 | 0.05 | 0.18 | -0.05 | 0.22 | 0.12
PCC $\uparrow$ | 0.75 | 0.54 | 0.46 | 0.24 | 0.63 | 0.36 | 0.42 | 0.25 | 0.48 | 0.24 | 0.47 | 0.37
Cliff | MSE $\downarrow$ | 4.36 | 2.88 | 5.94 | 2.90 | 5.34 | 2.98 | 5.96 | 3.07 | 5.51 | 3.17 | 4.47 | 3.06
MAE $\downarrow$ | 1.62 | 1.24 | 1.79 | 1.33 | 1.78 | 1.34 | 1.89 | 1.35 | 1.72 | 1.33 | 1.55 | 1.30
$R^{2}$ $\uparrow$ | 0.61 | 0.09 | 0.47 | 0.08 | 0.53 | 0.06 | 0.47 | 0.05 | 0.51 | 0.00 | 0.60 | 0.04
PCC $\uparrow$ | 0.79 | 0.34 | 0.69 | 0.30 | 0.73 | 0.27 | 0.69 | 0.29 | 0.74 | 0.29 | 0.78 | 0.27
Figure 5:  TransMA multi-modal feature fusion process based on molecule 3D
Transformer and molecule Mamba under the scaffold and cliff data splitting
methods in Hela and RAW 264.7 cell lines. (A)The UMAP plot of the feature
fusion process under scaffold data splitting in Hela cell line. (B)The UMAP
plot of the feature fusion process under cliff data splitting in Hela cell
line. (C)The UMAP plot of the feature fusion process under scaffold data
splitting in RAW 264.7 cell lines. (D)The UMAP plot of the feature fusion
process under scaffold data splitting in RAW 264.7 cell line.

###  3.4 Ablation Experiment

TransMA’s ability to extract molecular features comes from both molecule 3D
Transformer and molecule Mamba. To demonstrate the superiority of multi-modal
feature extraction for molecules, a ablation experiment is designed: comparing
the accuracy of predicting transfection efficiency between TransMA, molecule
3D Transformer, and molecule Mamba under scaffold and cliff data splitting
methods in Hela and RAW 264.7 cell lines. In the ablation experiment, TransMA,
molecule 3D Transformer, and molecule Mamba ensure consistency in training set
and testing set, as well as parameter settings. Table 2 presents the
performance of TransMA, molecule 3D Transformer, and molecule Mamba in
predicting transfection efficiency under two different data splitting methods
in Hela and RAW 264.7 cell lines. The results demonstrate that TransMA
outperforms both molecule 3D Transformer and molecule Mamba in predicting
transfection efficiency. The ablation experiment indicates that the fusion of
molecule 3D Transformer and molecule Mamba for extracting multi-dimensional
molecular structural features is a key factor in achieving superior
performance in the task of transfection efficiency prediction.

Table 2:  The results of predicting LNPs transfection efficiency by TransMA, molecule 3D Transformer, and molecule Mamba under scaffold and cliff data splitting methods in Hela and RAW 264.7 cell lines are presented. The best results are highlighted in bold.  Splitting | Scaffold | cliff
---|---|---
Method | TransMA | Molecule Transformer | Molecule Mamba | TransMA | Molecule Transformer | Molecule Mamba
Hela | RAW | Hela | RAW | Hela | RAW | Hela | RAW | Hela | RAW | Hela | RAW
MSE $\downarrow$ | 3.64 | 1.63 | 5.13 | 1.77 | 6.38 | 1.90 | 4.36 | 2.88 | 4.63 | 3.02 | 5.29 | 3.00
MAE $\downarrow$ | 1.57 | 0.98 | 1.88 | 1.02 | 2.17 | 1.07 | 1.62 | 1.24 | 1.62 | 1.29 | 1.74 | 1.30
$R^{2}$ $\uparrow$ | 0.49 | 0.23 | 0.56 | 0.17 | 0.11 | 0.10 | 0.61 | 0.09 | 0.59 | 0.05 | 0.53 | 0.06
PCC $\uparrow$ | 0.75 | 0.54 | 0.28 | 0.43 | 0.39 | 0.34 | 0.79 | 0.34 | 0.78 | 0.27 | 0.74 | 0.30

###  3.5 Analysis of model interpretability

In this section, the interpretability of the model is analyzed to verify
whether TransMA can identify the key atoms in ionizable lipids that affect
transfection efficiency. First, identifying the locations of key atoms in
ionizable lipids is fundamental for interpretability analysis. Therefore,
constructing transfection cliffs is employed to pinpoint key atoms.
Concurrently, the mol-attention mechanism block calculates the influence
scores of all atoms in ionizable lipids on transfection efficiency. If the
high-scoring atoms identified by the mol-attention mechanism align with the
key atoms found through the construction of transfection cliffs, it proves
that the model can accurately recognize the key atoms in ionizable lipids that
impact transfection efficiency.

####  3.5.1 Transfection cliffs

Transfection cliff molecular pairs refer to ionizable lipid molecules that are
structurally very similar but have significantly different transfection
efficiencies. The phenomenon of transfection cliffs exists in LNP datasets and
can be observed through UMAP plots where points that are close in space
exhibit large differences in mapped colors. Despite the structural similarity
of transfection cliff molecular pairs, there are differences in their atomic
compositions, which lead to significant differences in transfection
efficiency. Therefore, the atoms that differ between the structures of
transfection cliff molecular pairs are considered key atoms affecting
transfection efficiency. To quantify the similarity of ionizable lipid
molecules, similarity scores based on substructure similarity, scaffold
similarity, and SMILES string similarity are constructed[44]. Substructure
similarity and scaffold similarity are determined by calculating the Tanimoto
coefficient of the Extended Connectivity Fingerprints (ECFP)[45, 46] and the
Molecular ACCess System (MACCS) keys[47, 48] for the molecular graph
frameworks, respectively. SMILES string similarity is assessed by calculating
the Levenshtein distance[49] between the SMILES string representations of the
molecules. Equation 9 calculates the Tanimoto coefficient shown as:

$\text{Tanimoto coefficient}=\frac{c}{a+b-c},$ (9)

For substructure similarity, $a$ counts the number of set bits (1s) in
molecule A’s ECFP fingerprint, $b$ does the same for molecule B, and $c$
counts the number of common set bits in both fingerprints. For scaffold
similarity, $a$ counts set bits in molecule A’s MACCS keys, $b$ does the same
for molecule B, and $c$ counts common set bits in both molecules’ MACCS keys.
Equation 10 computes SMILES molecular similarity shown as:

$\text{SMILES string
similarity}=1-\frac{d(s_{1},s_{2})}{\max(|s_{1}|,|s_{2}|)},$ (10)

where $d(s1,s2)$ is the Levenshtein distance between the SMILES strings of two
molecules.

Similar molecule pairs are defined as pairs with a structural similarity
greater than 0.9. Equation 11 describes the structural similarity, which is
the average of substructure similarity, scaffold similarity, and SMILES string
similarity shown as:

$\text{structure similarity}=\frac{subs+scas+smis}{3},$ (11)

where $subs$ represents substructure similarity, $scas$ represents scaffold
similarity, and $smis$ represents SMILES string similarity. Since the
transfection efficiency in the LNPs dataset has been log2-transformed, similar
molecule pairs are identified as transfection cliff pairs when the
transfection difference exceeds 1. Equation 12 describes the transfection
difference shown as:

$transfection\ efficiencies=|\log_{10}(2^{{m_{2}-m1}})|,$ (12)

where $m_{1}$ and $m_{2}$ represent the transfection efficiencies of the two
molecules in the similar pair. A total of 4267 and 2104 transfection cliff
pairs are identified in the Hela and RAW 264.7 cell lines, respectively.
Figure 6 shows the scatter plots of structural similarity and transfection
difference for transfection cliff pairs in the Hela and RAW 264.7 cell lines.
From the figure, it can be seen that when the structural similarity of the
molecules is greater than 0.9, the transfection difference is significant,
reaching up to a thousand-fold or even ten thousand-fold difference.

Figure 6:  Scatter plots of structural similarity and transfection difference
for transfection cliff pairs in Hela and RAW 264.7 cell lines. (A) Scatter
plot of structural similarity and transfection difference for transfection
cliff pairs in the Hela cell line. (B) Scatter plot of structural similarity
and transfection difference for transfection cliff pairs in the RAW 264.7 cell
line.

####  3.5.2 Model interpretability

The transfection cliff phenomenon can help identify key atoms in ionizable
lipids that affect transfection efficiency. Figure 6 shows that even when the
molecular similarity of transfection cliff pairs is close to 1, the
transfection difference remains significant. This implies that the one or two
differing atoms between the two molecules in the transfection cliff pairs are
the cause of the large transfection difference. Therefore, the key atoms are
those differing atoms in the transfection cliff pairs. When predicting the
transfection efficiency of transfection cliff pairs, the mol-attention
mechanism block in TransMA calculates the attention scores for all atoms.
Figure 7 shows the attention scores calculated by the mol-attention mechanism
block and the key atoms of the transfection cliff pairs. It can be seen that
the atoms with the highest attention scores are precisely the key atoms
responsible for the differences in the transfection cliff pairs. For example,
in figure 7, the similarity of the transfection cliff pair is 0.91, and the
transfection difference is 1.69. This indicates that the transfection
efficiency of one molecule is $10^{1.69}$ times higher than the other. The two
molecules differ by only one carbon atom (C) and one nitrogen atom (N) in
structure. The attention scores corresponding to the differing C and N atoms
are 0.84 and 0.86, significantly higher than the attention scores for other
atoms. Therefore, TransMA can identify the key atoms affecting transfection
efficiency, demonstrating the interpretability of the model.

Figure 7:  The scatter plot of attention scores distribution and key atoms
display for transfection cliff pairs. The attention score ranges from 0 to 1.
The transfection difference value represents the multiplicative difference in
transfection efficiency between two molecules, expressed in powers of 10. For
instance, a transfection difference of 3.89 implies that the transfection
efficiency of one molecule is $10^{3.89}$ times higher than the other. The
structural similarity ranges from 0 to 1.

###  3.6 External test

While TransMA achieves accurate predictions of transfection efficiency on the
LNPs dataset, we aim to evaluate its generalization capability by testing it
on an external dataset without additional training. We compile an external
dataset based on published studies, consisting of 15 LNPs with ionizable
lipids as the only variable and their transfection efficiency levels. Among
these, 11 ionizable lipids which is from different combinations of 6 head
groups, 1 ester linkage, and 5 hydrophobic tails shown in figure 8 and their
transfection efficiencies[50] are obtained from Balb/c mice using a fixed LNP
formulation. The remaining 4 ionizable lipids and their transfection
efficiencies[51, 52, 53] are derived from different cell lines. Since the
TransMA model is also trained on different cell lines, it possesses the
capability to predict transfection efficiency across various cell lines.

Figure 8: Molecular structures of the 11 ionizable lipids

Figure 9 shows the results of TransMA directly testing on an external dataset
without training. The predicted maximum and minimum transfection efficiencies
of lipids 1-11 by TransMA are consistent with the true values, and the order
of predicted values is similar to that of true values. For lipids 12 and 13,
the prediction error is less than 1. As for lipids 14 and 15, the authors only
provided the transfection efficiency of lipid 14 being lower than that of
lipid 15 during collection, which is consistent with the prediction result.
TransMA achieves good predictions of transfection efficiency even without
training on the external dataset, demonstrating its strong generalization
capability.

Figure 9: Predicted results for the external dataset without training

##  4 Discussion and conclusion

In this work, we propose an explainable high-accuracy LNPs transfection
efficiency prediction model called TransMA. To achieve high-accuracy
predictions, TransMA employs a multi-modal molecular structure fusion
architecture, where the substructure molecule 3D Transformer extracts three-
dimensional spatial features of the molecule, and molecule mamba extracts one-
dimensional molecular features. To achieve model interpretability and feature
fusion, we design the mol-attention mechanism block. This block can integrate
multi-dimensional features and reveal atom-level structure-transfection
relationships based on the molecular channel attention mechanism. Compared
with advanced molecular graph convolutional networks and Transformer models,
TransMA achieves the highest accuracy in predicting transfection efficiency
under the scaffold and cliff data splitting methods in Hela and RAW 264.7 cell
lines. Moreover, we introduce transfection cliff pairs. The results
demonstrate that the atoms with high attention scores computed by the mol-
attention mechanism block correspond to the key atoms in the transfection
cliff pairs, showing that TransMA can identify the key atoms affecting
transfection efficiency. Additionally, the external test results indicate that
even on the untrained external dataset, TransMA’s predicted values maintain a
consistent order with the actual transfection efficiencies.

Although TransMA has demonstrated excellent performance in predicting
transfection efficiency, there are still limitations in prediction accuracy
due to the scarcity of LNPs datasets. Additionally, since ionizable lipids are
composed mainly of head and tail groups, different ionizable lipids tend to
have high structural similarity, making the occurrence of transfection cliffs
more likely. These transfection cliffs increase the difficulty for TransMA to
capture the mapping relationship between molecular structure and transfection
efficiency. Therefore, future work should focus on addressing the transfection
cliff phenomenon to further improve prediction accuracy.

## References

  * [1] Xucheng Hou, Tal Zaks, Robert Langer, and Yizhou Dong.  Lipid nanoparticles for mrna delivery.  Nature Reviews Materials, 6(12):1078–1094, 2021.
  * [2] Satoshi Uchida, Federico Perche, Chantal Pichon, and Horacio Cabral.  Nanomedicine-based approaches for mrna delivery.  Molecular Pharmaceutics, 17(10):3654–3684, 2020.
  * [3] Chaoyang Meng, Zhe Chen, Gang Li, Thomas Welte, and Haifa Shen.  Nanoplatforms for mrna therapeutics.  Advanced Therapeutics, 4(1):2000099, 2021.
  * [4] Makda S. Gebre, Luis A. Brito, Lisa H. Tostanoski, Darin K. Edwards, Andrea Carfi, and Dan H. Barouch.  Novel approaches for vaccine development.  Cell, 184(6):1589–1603, 2021.
  * [5] Yuhua Weng, Chunhui Li, Tongren Yang, Bo Hu, Mengjie Zhang, Shuai Guo, Haihua Xiao, Xing-Jie Liang, and Yuanyu Huang.  The challenge and prospect of mrna therapeutics landscape.  Biotechnology Advances, 40:107534, 2020.
  * [6] Jeonghwan Kim, Yulia Eygeris, Mohit Gupta, and Gaurav Sahay.  Self-assembled mrna vaccines.  Advanced Drug Delivery Reviews, 170:83–112, 2021.
  * [7] Wei Li, Yijia Li, Jingjiao Li, Junli Meng, Ziqiong Jiang, Chen Yang, Yixing Wen, Shuai Liu, Xingdi Cheng, Shiwei Mi, Yuanyuan zhao, Lei Miao, and Xueguang Lu.  All-trans-retinoic acid-adjuvanted mrna vaccine induces mucosal anti-tumor immune responses for treating colorectal cancer.  Advanced Science, n/a(n/a):2309770.
  * [8] Junjiao Song, Hongquan Li, Yanfang Liu, Xinrong Li, Qili Shi, Qun-Ying Lei, Weiguo Hu, Shenglin Huang, Zhiao Chen, and Xianghuo He.  Aldolase a accelerates cancer progression by modulating mrna translation and protein biosynthesis via noncanonical mechanisms.  Advanced Science, 10(26):2302425, 2023.
  * [9] Daisy Yi Ding, Yuhui Zhang, Yuan Jia, and Jiuzhi Sun.  Machine learning-guided lipid nanoparticle design for mrna delivery.  In ICML 2023 CompBio Workshop, 2023.
  * [10] Shugang Qin, Xiaoshan Tang, Yuting Chen, Kepan Chen, Na Fan, Wen Xiao, Qian Zheng, Guohong Li, Yuqing Teng, Min Wu, and Xiangrong Song.  mrna-based therapeutics: powerful and versatile tools to combat diseases.  Signal Transduction and Targeted Therapy, 7(1):166, 2022.
  * [11] Rumiana Tenchov, Robert Bird, Allison E. Curtze, and Qiongqiong Zhou.  Lipid nanoparticles-from liposomes to mrna vaccine delivery, a landscape of research diversity and advancement.  ACS Nano, 15(11):16982–17015, 2021.
  * [12] Cheng Huang, Xing Duan, Jichao Wang, Qingqing Tian, Yangmei Ren, Kepan Chen, Zongliang Zhang, Yuanyou Li, Yunyu Feng, Kunhong Zhong, Yuelong Wang, Liangxue Zhou, Gang Guo, Xiangrong Song, and Aiping Tong.  Lipid nanoparticle delivery system for mrna encoding b7h3-redirected bispecific antibody displays potent antitumor effects on malignant tumors.  Advanced Science, 10(3):2205532, 2023.
  * [13] Anna Patrignani, Nicolo’ Schicchi, Francesca Calcagnoli, Elena Falchetti, Nino Ciampani, Giulio Argalia, and Antoine Mariani.  Acute myocarditis following comirnaty vaccination in a healthy man with previous sars-cov-2 infection.  Radiology Case Reports, 16:3321 – 3325, 2021.
  * [14] Sahabi Kabir Sulaiman, Muhammad Sale Musa, Fatimah Isma’il Tsiga-Ahmed, Abdulwahab Kabir Sulaiman, and Abdulaziz Tijjani Bako.  A systematic review and meta-analysis of the global prevalence and determinants of covid-19 vaccine acceptance and uptake in people living with hiv.  Nature Human Behaviour, 2023.
  * [15] M. Kim, M. Jeong, S. Hur, Y. Cho, J. Park, H. Jung, Y. Seo, H. A. Woo, K. T. Nam, K. Lee, and H. Lee.  Engineered ionizable lipid nanoparticles for targeted delivery of rna therapeutics into different types of cells in the liver.  Science Advances, 7(9):eabf4398, 2021.
  * [16] Khalid A. Hajj and Kathryn A. Whitehead.  Tools for translation: non-viral materials for therapeutic mrna delivery.  Nature Reviews Materials, 2(10):17056, 2017.
  * [17] Piotr S. Kowalski, Arnab Rudra, Lei Miao, and Daniel G. Anderson.  Delivering the messenger: Advances in technologies for therapeutic mrna delivery.  Molecular Therapy, 27(4):710–728, 2019.
  * [18] Lei Miao, Linxian Li, Yuxuan Huang, Derfogail Delcassian, Jasdave Chahal, Jinsong Han, Yunhua Shi, Kaitlyn Sadtler, Wenting Gao, Jiaqi Lin, Joshua C. Doloff, Robert Langer, and Daniel G. Anderson.  Delivery of mrna vaccines with heterocyclic lipids increases anti-tumor efficacy by sting-mediated immune cell activation.  Nature Biotechnology, 37(10):1174–1185, 2019.
  * [19] Bowen Li, Rajith Singh Manan, Shun-Qing Liang, Akiva Gordon, Allen Jiang, Andrew Varley, Guangping Gao, Robert Langer, Wen Xue, and Daniel Anderson.  Combinatorial design of nanoparticles for pulmonary mrna delivery and genome editing.  Nature Biotechnology, 41(10):1410–1415, 2023.
  * [20] Bowen Li, Idris O. Raji, Akiva G. R. Gordon, Lizhuang Sun, Theresa M. Raimondo, Favour A. Oladimeji, Allen Y. Jiang, Andrew Varley, Robert S. Langer, and Daniel G. Anderson.  Accelerating ionizable lipid discovery for mrna delivery using machine learning and combinatorial chemistry.  Nature Materials, 2024.
  * [21] Xuexiang Han, Hanwen Zhang, Kamila Butowska, Kelsey L. Swingle, Mohamad-Gabriel Alameh, Drew Weissman, and Michael J. Mitchell.  An ionizable lipid toolbox for rna delivery.  Nature Communications, 12(1):7233, 2021.
  * [22] Ning Meng and Dirk Grimm.  Membrane-destabilizing ionizable phospholipids: Novel components for organ-selective mrna delivery and crispr–cas gene editing.  Signal Transduction and Targeted Therapy, 6(1):206, 2021.
  * [23] Xu Yue, Ma Shihao, Cui Haotian, Chen Jingan, Xu Shufen, Wang Kevin, Varley Andrew, Lu Rick Xing Ze, Wang Bo, and Li Bowen.  Agile platform: A deep learning-powered approach to accelerate lnp development for mrna delivery.  bioRxiv, page 2023.06.01.543345, 2023.
  * [24] Kun Wu, Xiulong Yang, Zixu Wang, Na Li, Jialu Zhang, and Lizhuang Liu.  Data-balanced transformer for accelerated ionizable lipid nanoparticles screening in mRNA delivery.  Briefings in Bioinformatics, 25(3):bbae186, 04 2024.
  * [25] Seyone Chithrananda, Gabriel Grand, and Bharath Ramsundar.  Chemberta: Large-scale self-supervised pretraining for molecular property prediction.  CoRR, abs/2010.09885, 2020.
  * [26] José Jiménez-Luna, Francesca Grisoni, and Gisbert Schneider.  Drug discovery with explainable artificial intelligence.  Nature Machine Intelligence, 2(10):573–584, 2020.
  * [27] Chuizheng Meng, Loc Trinh, Nan Xu, James Enouen, and Yan Liu.  Interpretability and fairness evaluation of deep learning models on mimic-iv dataset.  Scientific Reports, 12(1):7166, 2022.
  * [28] Tianyuan Liu, Junyang Huang, Delun Luo, Liping Ren, Lin Ning, Jian Huang, Hao Lin, and Yang Zhang.  Cm-sirpred: Predicting chemically modified sirna efficiency based on multi-view learning strategy.  International Journal of Biological Macromolecules, 264:130638, 2024\.
  * [29] Tianyuan Liu, Huiyuan Qiao, Zixu Wang, Xinyan Yang, Xianrun Pan, Yu Yang, Xiucai Ye, Tetsuya Sakurai, Hao Lin, and Yang Zhang.  Codlncscape provides a self-enriching framework for the systematic collection and exploration of coding lncrnas.  Advanced Science, 11(22):2400009, 2024.
  * [30] Albert Gu and Tri Dao.  Mamba: Linear-time sequence modeling with selective state spaces, 2023\.
  * [31] Jun Ma, Feifei Li, and Bo Wang.  U-mamba: Enhancing long-range dependency for biomedical image segmentation, 2024.
  * [32] Ziyang Wang and Chao Ma.  Weak-mamba-unet: Visual mamba makes cnn and vit work better for scribble-based medical image segmentation, 2024.
  * [33] Yue Liu, Yunjie Tian, Yuzhong Zhao, Hongtian Yu, Lingxi Xie, Yaowei Wang, Qixiang Ye, and Yunfan Liu.  Vmamba: Visual state space model, 2024.
  * [34] Sijie Zhao, Hao Chen, Xueliang Zhang, Pengfeng Xiao, Lei Bai, and Wanli Ouyang.  Rs-mamba for large remote sensing image dense prediction, 2024.
  * [35] Albert Gu, Karan Goel, and Christopher Re.  Efficiently modeling long sequences with structured state spaces.  In International Conference on Learning Representations, 2022.
  * [36] Sanghyun Woo, Jongchan Park, Joon-Young Lee, and In So Kweon.  Cbam: Convolutional block attention module.  In Vittorio Ferrari, Martial Hebert, Cristian Sminchisescu, and Yair Weiss, editors, Computer Vision – ECCV 2018, pages 3–19, Cham, 2018. Springer International Publishing.
  * [37] Jongchan Park, Sanghyun Woo, Joon-Young Lee, and In So Kweon.  BAM: bottleneck attention module.  In British Machine Vision Conference 2018, BMVC 2018, Newcastle, UK, September 3-6, 2018, page 147. BMVA Press, 2018.
  * [38] Jie Hu, Li Shen, and Gang Sun.  Squeeze-and-excitation networks.  In 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 7132–7141, 2018.
  * [39] Yuyang Wang, Jianren Wang, Zhonglin Cao, and Amir Barati Farimani.  Molecular contrastive learning of representations via graph neural networks.  Nature Machine Intelligence, 4(3):279–287, 2022.
  * [40] Xiaomin Fang, Lihang Liu, Jieqiong Lei, Donglong He, Shanzhuo Zhang, Jingbo Zhou, Fan Wang, Hua Wu, and Haifeng Wang.  Geometry-enhanced molecular representation learning for property prediction.  Nature Machine Intelligence, 4(2):127–134, 2022.
  * [41] Derek van Tilborg, Alisa Alenicheva, and Francesca Grisoni.  Exposing the limitations of molecular machine learning with activity cliffs.  Journal of Chemical Information and Modeling, 62(23):5938–5951, 2022.
  * [42] Kevin Yang, Kyle Swanson, Wengong Jin, Connor Coley, Philipp Eiden, Hua Gao, Angel Guzman-Perez, Timothy Hopper, Brian Kelley, Miriam Mathea, Andrew Palmer, Volker Settels, Tommi Jaakkola, Klavs Jensen, and Regina Barzilay.  Analyzing learned molecular representations for property prediction.  Journal of Chemical Information and Modeling, 59(8):3370–3388, 2019\.
  * [43] Gengmo Zhou, Zhifeng Gao, Qiankun Ding, Hang Zheng, Hongteng Xu, Zhewei Wei, Linfeng Zhang, and Guolin Ke.  Uni-mol: A universal 3d molecular representation learning framework.  In The Eleventh International Conference on Learning Representations, 2023.
  * [44] Derek van Tilborg, Alisa Alenicheva, and Francesca Grisoni.  Exposing the limitations of molecular machine learning with activity cliffs.  Journal of Chemical Information and Modeling, 62(23):5938–5951, 2022.
  * [45] David Rogers and Mathew Hahn.  Extended-connectivity fingerprints.  Journal of Chemical Information and Modeling, 50(5):742–754, 2010\.  PMID: 20426451.
  * [46] Ryosuke Asahara and Tomoyuki Miyao.  Extended connectivity fingerprints as a chemical reaction representation for enantioselective organophosphorus-catalyzed asymmetric reaction prediction.  ACS Omega, 7(30):26952–26964, 2022.
  * [47] Xueying Yuan, Lifeng Li, Zenan Shi, Hong Liang, Shuhua Li, and Zhiwei Qiao.  Molecular-fingerprint machine-learning-assisted design and prediction for high-performance mofs for capture of nmhcs from air.  Advanced Powder Materials, 1(3):100026, 2022.
  * [48] Umit V. Ucak, Taek Kang, Junsu Ko, and Juyong Lee.  Substructure-based neural machine translation for retrosynthetic prediction.  Journal of Cheminformatics, 13(1):4, 2021.
  * [49] Bonnie Berger, Michael S. Waterman, and Yun William Yu.  Levenshtein distance, sequence comparison and biological database search.  IEEE Transactions on Information Theory, 67(6):3287–3294, 2021\.
  * [50] Ziwei You.  Design and synthesis of ionizable cationic lipids and their application in personalized antitumor vaccine.  Master’s thesis, Central China Normal University, 2023.
  * [51] Shufen Xu, Yue Xu, Nicholas C. Solek, Jingan Chen, Fanglin Gong, Andrew James Varley, Alex Golubovic, Anni Pan, Songtao Dong, Gang Zheng, and Bowen Li.  Tumor-tailored ionizable lipid nanoparticles facilitate il-12 circular rna delivery for enhanced lung cancer immunotherapy.  Advanced Materials, n/a(n/a):2400307.
  * [52] Afsane Radmand, Melissa P. Lokugamage, Hyejin Kim, Curtis Dobrowolski, Ryan Zenhausern, David Loughrey, Sebastian G. Huayamares, Marine Z. C. Hatit, Huanzhen Ni, Ada Del Cid, Alejandro J. Da Silva Sanchez, Kalina Paunovska, Elisa Schrader Echeverri, Aram Shajii, Hannah Peck, Philip J. Santangelo, and James E. Dahlman.  The transcriptional response to lung-targeting lipid nanoparticles in vivo.  Nano Letters, 23(3):993–1002, 2023.  PMID: 36701517.
  * [53] Daisy Yi Ding, Yuhui Zhang, Yuan Jia, and Jiuzhi Sun.  Machine learning-guided lipid nanoparticle design for mrna delivery.  In ICML 2023 CompBio Workshop, 2023.


# Mobile 3D Printing Robot Simulation with Viscoelastic Fluids

Uljad Berdica1,2∗, Yuewei Fu1∗, Yuchen Liu1∗, Emmanouil Angelidis3,4∗, Chen
Feng1†
https://ai4ce.github.io/M3DP-Sim ∗ equal contribution.1New York University,
Brooklyn, NY 11201, USA {yf1236, yl5680<EMAIL_ADDRESS>York University
Abu Dhabi, UAE<EMAIL_ADDRESS>of Neuromorphic Computing, fortiss -
Research Institute of the Free State of Bavaria<EMAIL_ADDRESS>of
Robotics, Artificial Intelligence and Embedded Systems, Technical University
of Munich†Chen Feng is the corresponding author.

###### Abstract

The system design and algorithm development of mobile 3D printing robots need
a realistic simulation. They require a mobile robot simulation platform to
interoperate with a physics-based material simulation for handling
interactions between the time-variant deformable 3D printing materials and
other simulated rigid bodies in the environment, which is not available for
roboticists yet. To bridge this gap and enable the real-time simulation of
mobile 3D printing processes, we develop a simulation framework that includes
particle-based viscoelastic fluid simulation and particle-to-mesh conversion
in the widely adopted Gazebo robotics simulator, avoiding the bottlenecks of
traditional additive manufacturing simulation approaches. This framework is
the first of its kind that enables the simulation of robot arms or mobile
manipulators together with viscoelastic fluids. The method is tested using
various material properties and multiple collaborating robots to demonstrate
its simulation ability for the robots to plan and control the printhead
trajectories and to visually sense at the same time the printed fluid
materials as a free-form mesh. The scalability as a function of available
material particles in the simulation was also studied. A simulation with an
average of 5 FPS was achieved on a regular desktop computer.

##  I Introduction

Additive manufacturing (AM) is an increasingly popular production paradigm
that replaces labour intensive processes with highly integrated, numerically
controlled technologies that facilitate high adaptability to changing trends
and increasing complexity demands in product design [1, 2]. With the increased
use of such technology, several major issues arise like print size limitation,
print process duration [3] and different materials integration [4] along with
print quality [5]. Advances in robotic technologies [5, 6, 7] have catalyzed
the development of mobile 3D printing machinery which can overcome the
traditional limitations of gantry-based printing and build arbitrarily large
structures through coordinated collaborative processes between printer robots
[8, 9, 5].

AM projects of such scale require a higher degree of planning for pivotal
design decisions like required support [10], reliable material behavior that
maintains the printing patterns after extrusion [2] and error control for the
printer robots and printing materials [4]. Therefore, simulation tools for
mobile 3D printing are of paramount importance since they help address these
issues before production. Current work is either focused on the geometric
simulation of collaborative 3D printing [11, 9] or Finite Elements Method
(FEM) that balance the trade-offs between the costly meshing [12, 13, 14] and
signed distance field [15] computations.The simulations also lack real time
visualization of the layer by layer deposition and temporal change in material
behavior [12].

Figure 1: The viscoelastic fluid simulation with a robotic arm printing a
squared wall (left) and a Turtlebot3 printing a simple curved geometry
(right), each using different material properties in the same environment.

The contributions of this paper are as follows:

  * -

The proposal of a robotic AM simulation method that for the first time enables
simulated robots to visually perceive, plan, and control viscoelastic fluids,
bridging the gap between robotic and AM simulations.

  * -

The conduction of experiments demonstrating the effectiveness and near real-
time simulation frame rates of this framework (Figure 1), implemented as a
Gazebo plugin that will be open-sourced for research and education.

Note that this paper’s goal is to introduce this simulation method for
realistic mobile 3D printing robot simulations. We hope the community can take
advantage of such new simulation abilities to develop and test new mobile 3D
printing robot planning, control, and perception algorithms, which are also
the future work of this paper.

##  II Related Work

###  II-A Mobile 3D printing

The recent works in mobile 3D printing robotics have further established
mobile additive manufacturing as a valid solution to the limitations of static
AM [5, 11, 16, 3, 17]. Marques et al. [5] presented a mobile fused deposition
modelling (FDM) printer along with a circuit design for the remote control of
the robot through internet. The work by The Singapore Center for 3D printing
utilized a simulation featuring only two robots concurrently working on the
same project and was done on OpenRAVE without considering the time evolution
of material behavior.

McPherson and Zhou [9] propose a "chunk-by-chunk" slicer method to separate
the design into multiple parts that can be concurrently printed by a large
number of robots. Their simulator can only visualize the geometry of the
printing process; the extruded material is represented by a collection of
cylinders organized along the interpolated trajectory. This guided the work
presented in this paper toward a more granular material physics based
approach.

###  II-B AM Simulation and Analysis

Simulations use the predictability of AM processes to check the feasibility of
the print, and control errors that may arise during production [1, 10].

Finite Elements Analysis of designs was the most used approach both in
commercial products like Ansys and COMSOL Multiphysics® [18, 19, 20, 21] and
in academic works [22, 13, 14, 23]. Chen [12] pointed out that meshing
(an-/isotropic) is the computational bottleneck of FEA methods which require
prioritizing resolution for regions of interest.

The works of Biegler et al. [24] and Lee et al. [25] present a challenge that
arises from the accumulation of error when printing thin layers on top of one
another. Consequently, the physical libraries chosen for this work use a
particle-based approach to simulate the extrusion and deposition of the
material, bypassing the costly meshing of the entire design and taking into
account the change in fluid properties like surface tension and viscosity. A
marching cubes methodology similar to that of Cattenone et al. [23] was used
to mesh the particles, thus not relying on the computational expenses of
meshing to perform the real time simulation.

###  II-C Robotic simulation platform

The process of joining the high fidelity particle-based material simulation
with the robot simulation environment is unprecedented in the community and it
therefore requires a platform that ensures seamless integration and promotes
the use of our tool by the robotics and the AM community.

The comprehensive reviews of Nogueira [26] and Santos Pessoa de Melo et al.
[27] cemented the authors’ belief that the ease of interfacing with the Robot
Operating System (ROS) [28] constitutes a main evaluation criterion in the
choice between potential platforms like V-REP, ARGoS, Unity and OpenRave.
Surveys in the academic community also reveal that the ratio of participants
who use Gazebo as the main tool for simulation is twice larger than that of
other listed software in the survey [27].

On their introductory paper for Gazebo, Koenig and Howard [29] describe that
Gazebo can have significant computational costs when simulating a number of
robots exceeding the order of tens. Alternatively, V-REP has shown less CPU
usage than Gazebo [26] and it features more mesh editing capabilities. While
Gazebo exhibits these drawbacks, it was ultimately chosen due to its seamless
integration with ROS and the community and the open-source license. The
computations overhead is partly made up for by the external libraries used in
our simulation and the lack of mesh customization capabilities is accounted
for by a new external meshing methodology which offers more relevant options
to the context of the AM simulation.

TABLE I: The comparison summary of our method vs. related robotic and AM
simulation tools

|  Fluid simulation |  Material behavior |  Multi robot control |  Real time simulation |  Rigid body physics |  Mesh pipeline |
---|---|---|---|---|---|---|---
SOFA[30] | SPH | + |  | + | + |  |
SPlisHSPlasH[31] | SPH | + |  | + |  |  |
McPherson et al.[9] |  |  | + |  | + |  |
V-REP[32] |  | + | + | + | + |  |
Fluid Engine Dev[33] | SPH | + |  | + |  | + |
Hu et al. [34] | MPM | + |  | + | + |  |
AM in Ansys/COMSOL[18, 21] |  | + |  |  | + | + |
Cattenone et al.[23] |  | + |  |  |  | + |
Our tool | SPH | + | + | + | + | + |

The main Fluid Simulation method is listed when the tool has a fluid
simulation capability that is different from FEM of solid meshes. Real time
refers to whether the tool has the capacity to render AM processes in real
time. Material behavior refers to non-fluid material properties, e.g.,
elasticity.

###  II-D Material Simulation and Meshing Methods

Physics-based material simulation is at the core of the proposed methodology
in this paper. The cost efficiency and the relevancy of the methods is just as
important as the ability of the user to interact with all the layers of the
simulation starting from the solver and constraints to the visual
representation. The open-source Simulation Open Framework Architecture (SOFA)
project uses XML files that allow the user to interface with all the aspects
of the simulation before and during run-time [30]. A similar run-time
versatility approach was taken in the meshing process from the location of the
simulated particles.

The mathematical models of viscoelastic fluids by [35, 36] directed the
attention to Smooth Particle Hydrodynamics (SPH) as one of the most viable
solutions for a wide range of materials.SPlisHSPlasH [37, 38, 39, 40] was
chosen due to its flexibility in defining the physical properties of the
simulation (i.e. resolution, viscosity, stiffness etc) and choice of solvers.

Our SPlisHSPlasH plugin to Gazebo uses an XML file description that
encapsulates all these different properties. The method described in [41] was
used to sample the rigid boundaries with particles that interact with the
fluid particles based on the SPH approach. The forces applied on the boundary
particles are subsequently applied on the rigid bodies. The numerical
stability of the simulation is guaranteed by the CFL condition, which ensures
that the particle with the maximum velocity does not cover a distance longer
that twice the particle radius within a timestep, thus preventing particles
from overlapping and crossing rigid boundaries. It was found that timesteps at
around 1ms were enough to guarantee numerical stability both for the robotics
and viscoelastic simulation. The implicit DFSPH method [42] was employed for
the simulation. The highly efficient Moving Least Squares Material Point
Method (MPM) was also investigated [34, 43, 44, 45] as an alternative to SPH
methods.

Meshing is the final phase of this paper’s proposed method as it shows the
final simulated print result. [46, 47, 48, 49, 50] were very helpful in
designing a meshing pipeline that starts with generating a level set and
integrating it with a marching cubes algorithm to extract the surface of free
moving particles. Doyub Kim’s work on the Fluid Engine Development book [33]
provided a particles-to-object solution that allowed for the generation of an
.obj file from any given frame and import it to the simulation environment.
Thus, this new methodology bypasses the costly meshing bottlenecks pointed out
by Chen [12].

Figure 2: The workflow of our simulation approach. This describes the run from
the XML Gazebo world file to the SPH and Gazebo Wrapper code which is then
rendered in the Gazebo GUI (gzclient) and meshed in the final object.

##  III Method

Figure 2 shows the pipeline of the proposed simulation method. The Gazebo
scene file provides the initial conditions which are in turn fed to the SPH
library that calculates the new particle positions, as well as the forces
applied on the rigid bodies for each step. The particles position is used to
render the particles using an object-oriented rendering engine that creates
rigid entities of variable size.

###  III-A SPH Overview

Smoothed particle hydrodynamics (SPH) [35, 36] is at the core of this
viscoelastic material simulation. The main equations to be solved are the well
known equations of continuity for incompressible flow and the 3D Navier-Stokes
equation in Lagrangian coordinates as expressed by Koschier et al. [31].

In the SPH formalism, the fluid is discretized to an ensemble of particles
representing the physical space, each equipped with a Smoothing Kernel. Each
spatial quantity of the particle (f.e. mass, pressure) is computed as a
weighted sum of the quantities of its neighbouring particles, within a radius
of interaction. The closer the neighbouring particle, is located, the greater
the influence it has on the properties of the particle into consideration as
illustrated in Figure 3. More formally let f be any particle quantity (such as
mass, pressure) to be computed, $\left\langle f(x)\right\rangle$. This
function can be rewritten as a the convolution with the Dirac-$\delta$
function over a specified domain $\Omega$:

$\left\langle
f(x)\right\rangle=\int_{\Omega}f({x}^{\prime})\delta(x-{x}^{\prime})d{x}^{\prime}.$
(1)

Figure 3: Visualization of the averaging process for a property f using a
kernel W(r,h). Source: Adapted from [51].

In SPH, the identity (3) can be approximated by replacing the Dirac-$\delta$
with any kernel $W(\vec{r},h)$ as $h\rightarrow 0$. h is the variable radius
over which the physical properties are approximated. Equation (1) can be
discretized to compute any property of particle i described by a function
$f_{i}$ over a domain of a given number of particles.

For a thorough SPH review, we refer the reader to [31].

###  III-B Meshing

Figure 4: Pipeline of the meshing process from the particle data to the obj
file that can be put back in the running simulation. Marching Cubes
implementation by Kim [33].

The meshing process takes the Cartesian coordinates of the fluid model
calculated in SPH and creates an implicit surface level set whose boundaries
is determined by the points that satisfy an equation $F(x,y,z)$. Through the
marching cubes algorithm, a uniform-density polygonal mesh can be extracted
from the implicitly defined surface [52]. Once the object file is generated,
it is packaged in a Simulation Description Format (SDF) file which is then
inserted back to the simulation during run time. Figure 4 shows the flow of
the particle data from the simulation to the explicit triangulation functions
and back to the Gazebo environment.

As Zhu and Bridson [53] explained, a distance function defines the surface
that can generate an explicit mesh through the marching cubes algorithm:

$\phi(x)=\left|x-\overline{x}\right|-\overline{r},$ (2)

$\overline{r}$ ${}={}$ | weighted average of the particles radii,
---|---
$\overline{x}$ ${}={}$ | weighted average of the neighboring particles position.

In [53], the weight is assigned through:

$w_{i}=\frac{k(\left|x-x_{i}\right|/R)}{\sum_{j}k(\left|x-x_{j}\right|/R)},$
(3)

$R$ ${}={}$ | radius of the neighborhood around particle $x_{i}$,
---|---
$k$ ${}={}$ | kernel with tails that drop to zero.

The choice of kernel is non-trivial for special concave regions or small holes
since it can lead to the average of the particle being outside the surface
[53]. The simulation of the interaction of material with different densities
poses another challenge which can be resolved by the using an anisotropic
kernel to capture the densities more accurately at the expense of higher
computational costs [47]. The application of Laplacian smoothing and the
principal components analysis performed are covered in detail on [50, 54].
Other faster, less accurate algorithms like spherical and blobby do not have
the averaging step Zhu-Bridson and Anisotropic methods have. Through the
process shown in Figure 4, the user can choose between the methods based on
their simulation needs.

###  III-C Integration

The libraries for the fluid simulation and the meshing of the print results
are integrated in the workflow in Figure 2. The XML file reads the Gazebo
world environment configuration and the initial conditions for the fluid
simulation. Gazebo computes the positions of the rigid bodies and updates them
on the SPH simulator side. This way the robotic simulation is influenced by
the presence of the viscoelastic fluid particles. The Gazebo Messages library
helps communicate the Gazebo topics through Google protobuf messages in the
GazeboGenerator and MsgFactory class. The node connected to the simulation
method will keep sending messages including positions of all the particles.
The Subscriber node calls the rendering function through the Object-Oriented
Graphics Rendering Engine.

In Gazebo, a particle emitter is attached to a cylindrical model which allows
for a simplified representation of the print head attached to the robot. The
cylindrical model is controlled by a gazebo::physics::model plugin. The
movement of this model is instrumental in checking the printing patterns and
different levels of curvature and thickness. The rigid boundaries are sampled
by particles that behave like the viscoelastic material according to [41]
ensuring that interaction forces are applied on the boundary particles and on
the rigid bodies. This ensures that the collisions between the robot and the
viscoelastic material are handled properly.

##  IV Experiments

To demonstrate the validity of the proposed approach, a simulation is used to
show how the results are affected by different material properties.
Furthermore, the simulation plugin is integrated with various robots in order
to connect the viscoelastic material simulation to mobile 3D printing
simulations. A performance analysis is carried out to demonstrate the trade-
offs of the approach proposed in this paper which allows the user to
prioritize based on the objective of their simulation.

###  IV-A Simple particle dropping demonstration

As seen in the literature review, certain simulations may require little to no
realistic material behavior and are focused on the geometrical layout of the
layer by layer printing process. Therefore, the first experiment consists in
printing a cylindrical wall using only rigid entities under the normal gravity
force. Figure 5 shows the simulation of a simple block-dropping print with the
emitter head represented by a small cylinder.

Figure 5: The printing material is rendered into discrete rigid spheres. This
does not take into account the viscoelastic properties of the particle system.

###  IV-B Simulating Different Material Behavior

After having demonstrated the ability to show the layer by layer printing
using rigid bodies, this experiments shows the capability of the SPH libraries
to direct the rendered entities to behave as a fluid with time dependent
viscosity. Viscosity coefficient $\mu$, can be calculated from equation
(LABEL:NavStok) as product of kinematic viscosity coefficient $\eta$ and
density $\rho$ .

Figure 6: Effect of different surface tension and viscosity.

Figure 6 shows the capability to simulate a wide array of materials. In
addition to surface tension and viscosity, the SPH libraries account for
elasticity, vorticity and drag force.

###  IV-C Meshing of the simulated particles

The coordinates of the center of mass of each rendered entity (e.g. spheres)
are passed to the meshing function with a user-adjustable period. All the
emitted particles at a given step are used to generate the mesh which is then
fed back into the simulation environment through an internally generated SDF
file. Figure 7 shows the process for one layer.

Figure 7: With the given initial conditions, the emitter follows a
parameterized path and extrudes the particles that are meshed.

###  IV-D Integration with Mobile Robotics Simulation

The particle-emitting cylinder is replaced by a moving robot model by
attaching the emitter directly to the robot’s arm in its description file. A
Turtlebot Waffle Pi with openManipulator is used as the mobile printing
platform in the Gazebo environment (shown in Figure 1). OpenManipulator is
controlled using the movegroup interface of the Moveit! package. The Turtlebot
itself is controlled via odometry inputs to ROS. The Universal Robot Arm 10
(ur10) was also used to demonstrate the printing of the material by following
pre-determined waypoints.

Figure 8 left shows the printing process for the TurtleBot, where particles
are being continuously meshed during run time. Figure 8 right shows the result
of meshing a layer and starting the second while the robot mounted camera
perceives the print. This showcases the connection of the material simulation,
meshing pipeline and the mobile robot control for mobile 3D printing.

Figure 8: Left: mobile 3D printing TurtleBot following a pre-planned circular
trajectory. During printing, the emitted particles are meshed in given time
intervals. Particles and meshes can co-exist in the simulation. Right: robot
perception simulation of mobile 3D printing. The emitted particles are meshed
using spherical method during the printing. The printing result is observed by
a simulated camera on the mobile robot.

###  IV-E Multiple Robot Printing

In order for the mobile 3D printing simulation tool to be effective, the above
demonstration should be valid for multiple mobile robots. Figure 9 shows
multiple robots working collaboratively to print a structure. As demonstrated
in the experiments above, each robot is able to sense the other robots and the
printed material for collision avoidance between the robots themselves and the
printed mesh.

Figure 9: Multiple TurtleBots collaborating through the proposed framework.

###  IV-F Performance Analysis

Time taken per simulation step was measured over the number of particles
emitted using the ctime library. A baseline where particles are being
simulated in the back-end without any visualization is taken as a reference. A
significant increase in time per simulation step – from less than 100ms to
near 800ms – occurs when the total number of particles rendered as spheres
exceed 3000, resulting in an almost frozen simulation. The workaround for such
phenomenon is for the particles to be meshed, and then imported into the
gazebo GUI as .obj models. This method mitigates the steep increase in
simulation time, thus achieving stable simulation time around 200ms per step
while still running the layer-by-layer simulation as shown in Figure 10.

Figure 10: Performance chart comparing run time for different options.

##  V Conclusions

This work introduced a novel simulation approach for mobile 3D printing in
response to the most recent trends in mobile AM robotics and simulation. The
versatility of Smooth Particle Hydrodynamics for high viscosity materials is
combined with the the ubiquitous dynamic environment of Gazebo. This new
approach to simulation solves the computationally costly bottlenecks of mesh-
based simulation by adding an external meshing module which lets the user
choose the implicit surface generation method and the resolution. The mesh
along with the SPH computations are available as rigid entities in the Gazebo
multi-robot environment.

The validity of the proposed method is shown by connecting the mathematical
theory of SPH with the explicit meshing principles in a realistic physics
simulation model. The attachment of the simulation plugin to different robots
is shown along with their capability to sense and interact with these newly
generated entities in the scene.

Despite the promising results, the scalability of the Gazebo GUI for a number
of particles above the $4^{th}$ order of magnitude requires further work. The
most time-stable solution thus far has proven to be the meshing of the first
few thousand of particles with the meshing pipeline described in this paper
and then reusing the particles or simply using the point coordinates for
continuous meshing.

The authors hope that this approach will be used by researcher to help with
planning more aspects of their mobile 3D printing processes and inspire people
who lack access to resources to simulate and visualize a vast array of
printing materials in robotic additive manufacturing processes.

Future work could include the increase of the Gazebo GUI support for a larger
number of rigid entities than already demonstrated and the full implementation
of Material Point Methods which can replace the SPH libraries in the presented
workflow and extend the versatility of this tool.

## Acknowledgments

The research is supported by NSF CPS program under CMMI-1932187. Emmanouil
Angelidis is supported by the European Union’s Horizon 2020 research and
innovation programme under grant agreement No. 945539 (SGA3) Human Brain
Project.

## References

  * Paritala et al. [2017] P. K. Paritala, S. Manchikatla, and P. K. Yarlagadda, “Digital manufacturing-applications past, current, and future trends,” _Procedia engineering_ , vol. 174, pp. 982–991, 2017.
  * Rivera et al. [2020] R. G. Rivera, R. García Alvarado, A. Martínez-Rocamora, and F. Auat Cheein, “A comprehensive performance evaluation of different mobile manipulators used as displaceable 3d printers of building elements for the construction industry,” _Sustainability_ , vol. 12, no. 11, p. 4378, 2020\.
  * Currence et al. [2017] J. Currence, R. Morales-Ortega, J. Steck, and W. Zhou, “A floor power module for cooperative 3d printing,” _Proceedings of the 28th Annual International Solid Freeform Fabrication Symposium_ , 2017.
  * Oropallo and Piegl [2016] W. Oropallo and L. A. Piegl, “Ten challenges in 3d printing,” _Engineering with Computers_ , vol. 32, no. 1, pp. 135–148, 2016.
  * Marques et al. [2017] L. G. Marques, R. A. Williams, and W. Zhou, “A mobile 3d printer for cooperative 3d printing,” in _Proceeding of the 28th International Solid Freeform Fabrication Symposium_ , 2017, pp. 1645–1660.
  * Lloret et al. [2015] E. Lloret, A. R. Shahab, M. Linus, R. J. Flatt, F. Gramazio, M. Kohler, and S. Langenberg, “Complex concrete structures: Merging existing casting techniques with digital fabrication,” _Computer-Aided Design_ , vol. 60, pp. 40–49, 2015.
  * Bock [2015] T. Bock, “The future of construction automation: Technological disruption and the upcoming ubiquity of robotics,” _Automation in Construction_ , vol. 59, pp. 113–121, 2015.
  * Tiryaki et al. [2018] M. E. Tiryaki, X. Zhang, and Q.-C. Pham, “Printing-while-moving: a new paradigm for large-scale robotic 3d printing,” _arXiv preprint arXiv:1809.07940_ , 2018.
  * McPherson and Zhou [2018] J. McPherson and W. Zhou, “A chunk-based slicer for cooperative 3d printing,” _Rapid Prototyping Journal_ , vol. 24, 10 2018.
  * Ostrander et al. [2019] J. K. Ostrander, L. Ryan, S. Dhengre, C. McComb, T. W. Simpson, and N. A. Meisel, “A comparative study of virtual reality and computer-aided design to evaluate parts for additive manufacturing,” in _International Design Engineering Technical Conferences and Computers and Information in Engineering Conference_ , vol. 59186. American Society of Mechanical Engineers, 2019, p. V02AT03A029.
  * Zhang et al. [2018] X. Zhang, M. Li, J. Lim, Y. Weng, Y. W. D. Tay, H. Pham, and Q. C. Pham, “Large-scale 3d printing by a team of mobile robots,” _Automation in Construction_ , vol. 95, pp. 98–106, 08 2018.
  * Chen [2019] S. Chen, “Investigation of fem numerical simulation for the process of metal additive manufacturing in macro scale,” Ph.D. dissertation, Université de Lyon, 2019.
  * Alharbi et al. [2020] M. Alharbi, V. I. Patel _et al._ , “Simulation of uniaxial stress–strain response of 3d-printed polylactic acid by nonlinear finite element analysis,” _Applied Adhesion Science_ , vol. 8, no. 1, pp. 1–10, 2020.
  * Pal et al. [2014] D. Pal, N. Patil, K. Zeng, and B. Stucker, “An Integrated Approach to Additive Manufacturing Simulations Using Physics Based, Coupled Multiscale Process Modeling,” _Journal of Manufacturing Science and Engineering_ , vol. 136, no. 6, 10 2014, 061022. [Online]. Available: https://doi.org/10.1115/1.4028580
  * Barclay et al. [2016] J. Barclay, V. Dhokia, and A. Nassehi, “Additive manufacturing simulation using signed distance fields,” in _Sustainable Design and Manufacturing 2016_ , R. Setchi, R. J. Howlett, Y. Liu, and P. Theobald, Eds. Cham: Springer International Publishing, 2016, pp. 435–444.
  * [16] I. for Advanced Architecture of Catalonia. [Online]. Available: http://robots.iaac.net/
  * Steck et al. [2017] J. Steck, R. Morales-Ortega, J. Currence, and W. Zhou, “A mobile robot gripper for cooperative 3d printing,” _Proceedings of the 28th Annual International Solid Freeform Fabrication Symposium_ , 2017.
  * [18] “Additive manufacturing simulation.” [Online]. Available: https://www.ansys.com/products/structures/additive-manufacturing
  * [19] “Buy services: Fortissimo marketplace.” [Online]. Available: https://www.fortissimo-project.eu/buy-services
  * alp [2020] “Additive manufacturing,” May 2020. [Online]. Available: https://www.alphastarcorp.com/additive-manufacturing/
  * [21] “Simulating manufacturing processes with comsol multiphysics®.” [Online]. Available: https://www.comsol.com/video/simulating-manufacturing-processes-with-comsol-multiphysics-nov-29-2018
  * Zhang et al. [2020] Y. Zhang, F. Li, and D. Jia, “Residual stress and deformation analysis of lattice compressor impeller based on 3d printing simulation,” _Mechanics of Advanced Materials and Structures_ , pp. 1–15, 2020.
  * Cattenone et al. [2019] A. Cattenone, S. Morganti, G. Alaimo, and F. Auricchio, “Finite element analysis of additive manufacturing based on fused deposition modeling: distortions prediction and comparison with experimental data,” _Journal of Manufacturing Science and Engineering_ , vol. 141, no. 1, 2019.
  * Biegler et al. [2018] M. Biegler, B. Graf, and M. Rethmeier, “Assessing the predictive capability of numerical additive manufacturing simulations via in-situ distortion measurements on a lmd component during build-up,” _Procedia CIRP_ , vol. 74, pp. 158 – 162, 2018, 10th CIRP Conference on Photonic Technologies [LANE 2018]. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S2212827118308539
  * Lee et al. [2019] M.-H. Lee, S.-I. Chen, W.-H. Chen, and Y. Mao, “A layer-based mesh generator and scheme for 3d printing simulation,” _Computer Modeling in Engineering & Sciences_, vol. 120, no. 2, pp. 363–374, 2019.
  * Nogueira [2014] L. Nogueira, “Comparative analysis between gazebo and v-rep robotic simulators,” _Seminario Interno de Cognicao Artificial-SICA_ , vol. 2014, no. 5, 2014.
  * Santos Pessoa de Melo et al. [2019] M. Santos Pessoa de Melo, J. Gomes da Silva Neto, P. Jorge Lima da Silva, J. M. X. Natario Teixeira, and V. Teichrieb, “Analysis and comparison of robotics 3d simulators,” in _2019 21st Symposium on Virtual and Augmented Reality (SVR)_ , 2019, pp. 242–251.
  * Quigley et al. [2009] M. Quigley, K. Conley, B. Gerkey, J. Faust, T. Foote, J. Leibs, R. Wheeler, and A. Ng, “Ros: an open-source robot operating system,” vol. 3, 01 2009.
  * Koenig and Howard [2004] N. Koenig and A. Howard, “Design and use paradigms for gazebo, an open-source multi-robot simulator,” in _2004 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS) (IEEE Cat. No.04CH37566)_ , vol. 3, 2004, pp. 2149–2154 vol.3.
  * Faure et al. [2012] F. Faure, C. Duriez, H. Delingette, J. Allard, B. Gilles, S. Marchesseau, H. Talbot, H. Courtecuisse, G. Bousquet, I. Peterlik _et al._ , “Sofa: A multi-model framework for interactive physical simulation,” in _Soft tissue biomechanical modeling for computer assisted surgery_. Springer, 2012, pp. 283–321.
  * Koschier et al. [2019] D. Koschier, J. Bender, B. Solenthaler, and M. Teschner, “Smoothed Particle Hydrodynamics Techniques for the Physics Based Simulation of Fluids and Solids,” in _Eurographics 2019 - Tutorials_ , W. Jakob and E. Puppo, Eds. The Eurographics Association, 2019, pp. 1–41.
  * Freese et al. [2010] M. Freese, S. Singh, F. Ozaki, and N. Matsuhira, “Virtual robot experimentation platform v-rep: A versatile 3d robot simulator,” vol. 6472, 11 2010, pp. 51–62.
  * Kim [2017] D. Kim, _Fluid engine development_. CRC Press, 2017.
  * Hu et al. [2018] Y. Hu, Y. Fang, Z. Ge, Z. Qu, Y. Zhu, A. Pradhana, and C. Jiang, “A moving least squares material point method with displacement discontinuity and two-way rigid body coupling,” _ACM Transactions on Graphics (TOG)_ , vol. 37, no. 4, pp. 1–14, 2018.
  * Chang et al. [2009] Y. Chang, K. Bao, Y. Liu, J. Zhu, and E. Wu, “A particle-based method for viscoelastic fluids animation,” in _Proceedings of the 16th ACM Symposium on Virtual Reality Software and Technology_ , 2009, pp. 111–117.
  * Gerszewski et al. [2009] D. Gerszewski, H. Bhattacharya, and A. W. Bargteil, “A point-based method for animating elastoplastic solids,” in _Proceedings of the 2009 ACM SIGGRAPH/Eurographics Symposium on Computer Animation_ , 2009, pp. 133–138.
  * Bender et al. [2017] J. Bender, D. Koschier, T. Kugelstadt, and M. Weiler, “A micropolar material model for turbulent sph fluids,” in _Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Computer Animation_ , 2017, pp. 1–8.
  * Bender et al. [2019] J. Bender, T. Kugelstadt, M. Weiler, and D. Koschier, “Volume maps: An implicit boundary representation for sph,” in _Motion, Interaction and Games_ , 2019, pp. 1–10.
  * Weiler et al. [2016] M. Weiler, D. Koschier, and J. Bender, “Projective fluids,” in _Proceedings of the 9th International Conference on Motion in Games_ , 2016, pp. 79–84.
  * Weiler et al. [2018] M. Weiler, D. Koschier, M. Brand, and J. Bender, “A physically consistent implicit viscosity solver for sph fluids,” in _Computer Graphics Forum_ , vol. 37, no. 2. Wiley Online Library, 2018, pp. 145–155.
  * Akinci et al. [2012] N. Akinci, M. Ihmsen, G. Akinci, B. Solenthaler, and M. Teschner, “Versatile rigid-fluid coupling for incompressible sph,” _ACM Trans. Graph._ , vol. 31, no. 4, Jul. 2012. [Online]. Available: https://doi.org/10.1145/2185520.2185558
  * Bender and Koschier [2015] J. Bender and D. Koschier, “Divergence-free smoothed particle hydrodynamics,” in _Proceedings of the 2015 ACM SIGGRAPH/Eurographics Symposium on Computer Animation_. ACM, 2015.
  * Hu and Fang [2017] Y. Hu and Y. Fang, “An asynchronous material point method,” in _ACM SIGGRAPH 2017 Posters_ , 2017, pp. 1–2.
  * Stomakhin et al. [2013] A. Stomakhin, C. Schroeder, L. Chai, J. Teran, and A. Selle, “A material point method for snow simulation,” _ACM Transactions on Graphics (TOG)_ , vol. 32, no. 4, pp. 1–10, 2013.
  * Ram et al. [2015] D. Ram, T. Gast, C. Jiang, C. Schroeder, A. Stomakhin, J. Teran, and P. Kavehpour, “A material point method for viscoelastic fluids, foams and sponges,” in _Proceedings of the 14th ACM SIGGRAPH/Eurographics Symposium on Computer Animation_ , 2015, pp. 157–163.
  * Shi et al. [2019] W. Shi, P. X. Liu, and M. Zheng, “A mixed-depth visual rendering method for bleeding simulation,” _IEEE/CAA Journal of Automatica Sinica_ , vol. 6, no. 4, pp. 917–925, 2019.
  * Wang et al. [2017] X. Wang, X. Ban, Y. Zhang, Z. Pan, and S. Liu, “Anisotropic surface reconstruction for multiphase fluids,” in _2017 International Conference on Cyberworlds (CW)_. IEEE, 2017, pp. 118–125.
  * Yu et al. [2012] J. Yu, C. Wojtan, G. Turk, and C. Yap, “Explicit mesh surfaces for particle based fluids,” in _Computer Graphics Forum_ , vol. 31, no. 2pt4. Wiley Online Library, 2012, pp. 815–824.
  * Orthmann et al. [2010] J. Orthmann, M. Keller, and A. Kolb, “Topology-caching for dynamic particle volume raycasting.” in _VMV_ , 2010, pp. 147–154.
  * Wang et al. [2012] W. Wang, Z. Jiang, H. Qiu, and W. Li, “Real-time simulation of fluid scenes by smoothed particle hydrodynamics and marching cubes,” _Mathematical Problems in Engineering_ , vol. 2012, 2012.
  * Dai et al. [2016] Z. Dai, F. Wang, Y. Huang, K. Song, and A. Iio, “Sph-based numerical modeling for the post-failure behavior of the landslides triggered by the 2016 kumamoto earthquake,” _Geoenvironmental Disasters_ , vol. 3, no. 1, pp. 1–14, 2016.
  * Lorensen and Cline [1987] W. E. Lorensen and H. E. Cline, “Marching cubes: A high resolution 3d surface construction algorithm,” _ACM siggraph computer graphics_ , vol. 21, no. 4, pp. 163–169, 1987.
  * Zhu and Bridson [2005] Y. Zhu and R. Bridson, “Animating sand as a fluid,” _ACM Trans. Graph._ , vol. 24, no. 3, p. 965–972, Jul. 2005. [Online]. Available: https://doi.org/10.1145/1073204.1073298
  * Yu and Turk [2013] J. Yu and G. Turk, “Reconstructing surfaces of particle-based fluids using anisotropic kernels,” _ACM Transactions on Graphics (TOG)_ , vol. 32, no. 1, pp. 1–12, 2013.