# SimLOB: Learning Representations of Limited Order Book for Financial Market
Simulation

Yuanzhe Li
Department of Statistics and Data Science
Southern University of Science and Technology
<EMAIL_ADDRESS>
&Yue Wu11footnotemark: 1
Department of Statistics and Data Science
Southern University of Science and Technology
<EMAIL_ADDRESS>
&Peng Yang
Guangdong Provincial Key Laboratory of Brain-Inspired Intelligent Computation
Department of Statistics and Data Science
Department of Computer Science and Engineering
Southern University of Science and Technology
<EMAIL_ADDRESS>
Contributed equally to this work. The authors are ordered
alphabetically.Corresponding author.

###### Abstract

Financial market simulation (FMS) serves as a promising tool for understanding
market anomalies and the underlying trading behaviors. To ensure high-fidelity
simulations, it is crucial to calibrate the FMS model for generating data
closely resembling the observed market data. Previous efforts primarily
focused on calibrating the mid-price data, leading to essential information
loss of the market activities and thus biasing the calibrated model. The Limit
Order Book (LOB) data is the fundamental data fully capturing the market
micro-structure and is adopted by worldwide exchanges. However, LOB is not
applicable to existing calibration objective functions due to its tabular
structure not suitable for the vectorized input requirement. This paper
proposes to explicitly learn the vectorized representations of LOB with a
Transformer-based autoencoder. Then the latent vector, which captures the
major information of LOB, can be applied for calibration. Extensive
experiments show that the learned latent representation not only preserves the
non-linear auto-correlation in the temporal axis, but the precedence between
successive price levels of LOB. Besides, it is verified that the performance
of the representation learning stage is consistent with the downstream
calibration tasks. Thus, this work also progresses the FMS on LOB data, for
the first time.

##  1 Introduction

One important and challenging issue in financial market data is to discover
the causes of the abnormal market phenomena, e.g., flash crash [1, 2] and
spoofing [3]. Traditional machine learning methods are well-equipped to detect
the anomalies but lack of the financial interpretability of the causes [4].
Given an observed market time series data, financial market simulation (FMS)
aims to directly approximate its underlying data generative process, informed
by the trading rules of the real market [5]. By simulating the underlying
trading activities behind the observed market data, FMS is able to both
discover how the structure of the financial market changed over time and
explain the changes on the micro level [6, 7].

Generally, FMS models the market as a parameterized multi-agent system
$M(\mathbf{w})$ that mimics various types of traders and the exchange as
different agents. Each trading agent only interacts with the exchange agent
who implements the real trading rules of the market. Since 1990s, diverse
trading agents have been designed to simulate the chartists, fundamentalists,
momentum traders, high-frequency traders and so on [8], covering most of the
trading behaviors recognized from the real markets. By running the system,
those trading agents continuously submit orders to the exchange agent and the
exchange produces the simulated market data by matchmaking those orders.
Normally, the exchange will publicize the newest market data
$\hat{\mathbf{x}}(t)$ at a certain frequency (e.g., 1 second). Hence, the FMS
can be viewed as a market data generative process, denoted as
$M(\mathbf{w})=\mathbf{X}^{\mathbf{w}}_{T}=\\{\mathbf{x}^{\mathbf{w}}(t)\\}_{t=1}^{T}$
with any length of $T\in\mathbb{N}^{+}$.

To simulate any given observed financial market data
$\hat{\mathbf{X}}_{T}=\\{\hat{\mathbf{x}}(t)\\}_{t=1}^{T}$, FMS requires the
simulated data $\mathbf{X}^{\mathbf{w}}_{T}$ to closely resemble the observed
data $\hat{\mathbf{X}}_{T}$, which forces $M(\mathbf{w})$ approximating the
underlying data generative process of $\hat{\mathbf{X}}_{T}$. This largely
relies on the careful calibration of $M(\mathbf{w})$ by tuning the parameters
$\mathbf{w}$ to minimize the discrepancy between the two data sequences,
denoted as
$D(\hat{\mathbf{X}}_{T},M(\mathbf{w}))=D(\hat{\mathbf{X}}_{T},\mathbf{X}^{\mathbf{w}}_{T})$.
Unfortunately, the calibration problem is non-trivial due to the highly non-
linear interactions among the simulated agents [9]. In recent years, various
advanced methods have emerged from both the fields of optimization and
statistical inference [10]. The former minimizes the discrepancy using black-
box optimization methods [11], while the latter estimates the likelihood or
posterior of $\hat{\mathbf{X}}_{T}$ and $\mathbf{w}$ on the selected samples
where $D(\hat{\mathbf{X}}_{T},M(\mathbf{w}))<\epsilon$ [12].

In both ways, the calibration is mostly guided by the discrepancy.
Intuitively, the more information observed from the market is included in
$\hat{\mathbf{X}}_{T}$, the closer the approximation of $M(\mathbf{w})$ to the
underlying data generative process can be expected. In the literature, almost
all FMS works merely consider $\hat{\mathbf{X}}_{T}$ as the mid-price data
[10], ignoring the other important observable data information, e.g., trading
volume, bid/ask directions, and orders inter-arrival time. Consequently, the
calibrated $M(\mathbf{w})$ cannot capture the full dynamics of the real data
generative process where $\hat{\mathbf{X}}_{T}$ is generated.

To our best knowledge, this is the first work of calibrating FMS with respect
to the Limited Order Book (LOB) data, the fundamental market data widely
adopted in most of the world-class securities exchanges [13]. LOB is a complex
data structure for continuously recording all the untraded orders submitted
from all traders. In LOB, all the untraded orders fall into either ask side
(for selling) or bid side (for buying) according to their order directions. At
each side, the untraded orders are organized in the "price-first-time-second"
manner [14]. When a new order comes, the LOB is updated immediately in two
possible ways: if the new order’s price matches any opposite price of LOB, it
will be traded and the corresponding untraded orders will be deleted from LOB;
Otherwise, the new order will be untraded and inserted in its own side of LOB.
Note that, the untraded orders implicitly reflect the trading intentions of
the whole market. Besides, the immediate updates of LOB represent the mostly
fine-grained dynamics of the market micro-structure. On this basis, simulating
LOB data conceptually defines the closest approximation to the real data
generative process of the market scenario of interest.

Unfortunately, it is not straightforward to apply LOB into existing
discrepancy measures. Representative discrepancies like Euclidean distances
[15, 16] and probabilistic distances [11] all require $\hat{\mathbf{X}}_{T}$
and $\mathbf{X}^{\mathbf{w}}_{T}$ to be vectorized inputs, while the $T$ time
steps LOB is publicized in the form of normally a $10\times 4\times T$ matrix
data. That is, at each time step $t$, the best $10$ price levels on both sides
of LOB as well as the associated total volume of the untraded orders are
output as $\hat{\mathbf{x}}(t)\in\mathbb{R}^{10\times 4}$. More specifically,
$\hat{\mathbf{x}}(t)=[[p^{b}_{1}(t),v^{b}_{1}(t),p^{a}_{1}(t),v^{a}_{1}(t)];...;[p^{b}_{10}(t),v^{b}_{10}(t),p^{a}_{10}(t),v^{a}_{10}(t)]]$,
where $p_{i}$ and $v_{i}$ are for the $i$-th price level and the associated
total volume, $a$ and $b$ indicates the ask side and bid side (see Figure 1),
respectively. Existing FMS works that calibrate $1\times T$ mid-price can be
viewed as simplifying LOB by averaging the best ask price and the best bid
price at $t$-th step, i.e.,
$\hat{mp}(t)=\frac{p^{a}_{1}(t)-p^{b}_{1}(t)}{2}\in\mathbb{R}^{1}$, which
naturally lose much important information of the underlying data generative
process.

Figure 1: An illustration of the LOB data structure.

This paper aims to learn a low-dimensional vectorized representation of the
LOB data with the autoencoder, so that the well-established discrepancy
measures can still be easily adopted for calibration. Ideally, if the output
of the decoder closely resembles the input of the encoder, the latent vector
is believed to be an effective representation of the input LOB. The key
challenge lies in that the LOB involves not only the non-linear auto-
correlation in the temporal axis but also the precedence between successive
price levels at each time step. The literature witnesses increasing efforts on
analyzing LOB data with neural networks, which, however, neither adopt the
autoencoder framework nor explicitly study the vectorized representations. To
accurately represent these two properties in the latent vector, the proposed
encoder as well as the symmetric decoder contains three blocks, i.e., a fully
connected network for extracting features of the price precedence, a
Transformer stack for capturing the temporal auto-correlation, and another
fully connected network for dimension reduction.

Three groups of experiments have been conducted. The empirical findings are as
follows:

Finding 1: Calibrating LOB leads to better FMS than traditional mid-price
calibration;

Finding 2: The vectorized representations of various LOB can be effectively
learned, while the convolution layers adopted in existing works is not
effective;

Finding 3: Better representation of LOB, better calibration of FMS.

The rest of this paper is organized as follows. Section 2 introduces the
related works of FMS calibration and deep learning for LOB. Section 3
describes the proposed Transformer-based autoencoder in details. Section 4
reports 3 key empirical findings. The conclusion is drawn in Section 5.

##  2 Related Work

Calibration of FMS. Traditional FMS works focus on designing rule-based agents
to imitate various trading behaviors. Though several phenomena generally
existed in different stocks have been reproduced and explained by FMSs [17,
18], they received increasing criticisms for not being able to simulate any
specific time series data but only macro properties of the market [19]. Since
2015, several calibration objectives have been proposed for FMS [10]. A
straightforward way [15] is to simply calculate the mean square error (MSE)
between the simulated mid-price vector and the observed mid-price vector.
Methods of Simulated Moments calculate the statistical moments of both the
observed and simulated mid-price, and measure the distances between the two
vectors of moments as the discrepancy [16]. To better utilize the temporal
properties of the mid-price vector, some information-theoretic criteria based
discrepancies are proposed based on various time windows [20]. Francesco [21]
tries to compares the distribution distance between the simulated and observed
mid-price data. The Kolmogorov-Smirnov (K-S) test is generalized to multi-
variate data, while their empirical verification was only conducted on the
mid-price vector and the traded volume vector [11]. To summarize, existing
calibration functions all require the vectorized input format and related
works mainly calibrate to the mid-price data.

Deep Learning for LOB. Most of the studied tasks are to predict the price
movement [22], a binary classification problem. Since classification with
neural networks is known to be effective, the end-to-end solutions are
naturally considered for predicting the price movement from LOB and the
vectorized representation learning is not explicitly defined. In the contrary,
since the FMS models are mostly rule-based (especially the matchmaking rule of
updating the LOB), they are non-differentiable and cannot be trained end-to-
end with the feature extraction layers. Thus, the vectorized representation
learning need especial treatments in FMS.

Recent works propose to generate LOB data with deep learning, which is quite
similar to the goal of FMS [23]. However, they do not involve any validation
or calibration steps to force the generated data resembling any given observed
LOB, but only limited to follow some macro properties of the market [24, 25].
Some of them also are not informed by the real-world matchmaking rules for
updating LOB [26, 27]. These issues impose significant restrictions on the use
of these models as simulators, since one can hardly intervene the order
streams to do "what-if" test and reveal the micro-structured causes of certain
financial events [5, 28]. Furthermore, they explicitly require the ground
truth order streams as input to train the networks, which is not available in
the setting of FMS calibration problems.

Among them, the convolution layers and Long-Short Term Memory (LSTM) are the
mostly used architecture. Recently, the Transformer model is adopted instead
of LSTM, but the convolution layers are still kept [29, 30]. This work argues
for abandoning the convolution layer as it is empirically found ineffective to
learning from LOB.

##  3 Method

To feed the LOB data to the well-established calibration functions $D$ who
only accept the vectorized input, we propose to learn a function $f$ that can
represent the time series LOB data in a vectorized latent space. The
requirement of $f$ is two-fold. First, $f$ naturally constitutes a dimension
reduction process. Suppose $f$ receives a $\tau$ step LOB
$\hat{\mathbf{X}}_{\tau}\in\mathbb{R}^{10\times 4\times\tau}$, it is expected
that $f(\hat{\mathbf{X}}_{\tau})\in\mathbb{R}^{1\times\tilde{\tau}}$ with
$\tilde{\tau}\in\mathbb{N}^{+}$ and $\tilde{\tau}\ll 40\tau$. Second, within
the latent space of $\mathbb{R}^{1\times\tilde{\tau}}$, the information of LOB
can be largely preserved for downstream calibration tasks. One verification is
the existence of an inverse function
$g:\mathbb{R}^{1\times\tilde{\tau}}\rightarrow\mathbb{R}^{10\times
4\times\tau}$ such that $g(f(\hat{\mathbf{X}}_{\tau}))$ resembles the original
$\hat{\mathbf{X}}_{\tau}$. To this end, the autoencoder architecture is
considered, where $f$ constitutes the encoder network and $g$ is modeled by
the decoder. The latent vector $\mathbf{Z}\in\mathbb{R}^{1\times\tilde{\tau}}$
between the encoder and the decoder is the learned representation of
$\hat{\mathbf{X}}_{\tau}$, i.e., $\mathbf{Z}=f(\hat{\mathbf{X}}_{\tau})$, and
$g(\mathbf{Z})$ is the reconstructed data.

The Properties of LOB. The difficulty of the representation learning lies in
the specific properties of LOB. First, the time series LOB is essentially
auto-correlated between successive time steps, since $\hat{\mathbf{x}}_{t+1}$
is generated by updating $\hat{\mathbf{x}}(t)$ with incoming order streams
between the time interval $[t,t+1]$. Second, $\hat{\mathbf{x}}(t)$ comprises
10 levels of price and volume on both ask and bid sides. The price levels
should strictly follow the precedence that prices at lower levels are worse
than those at higher levels. That is, the bid/ask price at the $i$-th level
should be certainly larger/smaller than the $i+1$-th level, where
$i\in\mathbb{N}^{+}$ and $i\leq 9$. Hence, the encoder network should be able
to effectively handle both the non-linear temporal auto-correlation and the
precedence of price levels. On this basis, the proposed encoder is comprised
of three components, i.e., the feature extraction block as a fully connected
network (FCN, denoted as $\mathrm{FCN}_{1}$), a Transformer stack, and the
dimension reduction block as another FCN (denoted as $\mathrm{FCN}_{2}$). The
architecture of the encoder is depicted in Figure 2.

Figure 2: The architecture of the proposed encoder.

The Feature Extraction Block. The $\mathrm{FCN}_{1}$ is designed for
extracting useful features from LOB. Previous works that utilize machine
learning methods for processing LOB data often pre-process a set of hand-
crafted features to describe market dynamics [31]. Those works typically
consider the price levels as independent features and seldom deal with the
precedence relationship. Recent pioneer work [32] uses convolution neural
networks (CNNs) to automatically extract the features from LOB. The intuition
is that CNN may capture the precedence between the price levels through
convolutions like what has been done to the pixels. However, the 10 levels of
bid prices, bid volumes, ask prices, and ask volumes, have quite different
scales and meanings, which is non-trivial to be convolved adequately. This
work utilizes the FCN for feature extraction by
$\mathbf{h}_{0}=\mathrm{FCN}_{1}(\hat{\mathbf{X}}_{\tau})$. First, the LOB at
each time step is flattened from
$\hat{\mathbf{x}}(t)=[[p^{b}_{1}(t),v^{b}_{1}(t),p^{a}_{1}(t),v^{a}_{1}(t)];...;[p^{b}_{10}(t),v^{b}_{10}(t),p^{a}_{10}(t),v^{a}_{10}(t)]]$
to
$[p^{b}_{1}(t),v^{b}_{1}(t),p^{a}_{1}(t),v^{a}_{1}(t),...,p^{b}_{10}(t),v^{b}_{10}(t),p^{a}_{10}(t),v^{a}_{10}(t)]^{\mathsf{T}}$
using 1 linear layer. That is, the flattened $\hat{\mathbf{x}}(t)$ is into
$\mathbb{R}^{40\times 1}$. The times series LOB $\hat{\mathbf{X}}_{\tau}$ is
accordingly flattened in $\mathbb{R}^{40\times\tau}$ and then projected into
higher-dimensional space of $\mathbf{h}_{0}\in\mathbb{R}^{256\times\tau}$ to
obtain a more diverse range of features, using also 1 linear layer.

The Transformer Block. The second component of the encoder is the Transformer
stack with $L$ linked vanilla Transformers [33]. It aims to learn the temporal
auto-correlation between successive time steps with the multi-headed self-
attention (MSA), feed-forward layer (or say FCN), and layernorm (LN).
Specifically, at each $l$-th stacked Transformer, it computes

$\displaystyle\mathbf{h}_{l}^{\prime}$
$\displaystyle=\mathrm{MSA}(\mathrm{LN}(\mathbf{h}_{l-1}))+\mathbf{h}_{l-1},\quad$
$\displaystyle l=1,\ldots,L$ (1) $\displaystyle\mathbf{h}_{l}$
$\displaystyle=\mathrm{FCN}(\mathrm{LN}(\mathbf{h}_{l}^{\prime}))+\mathbf{h}_{l}^{\prime},\quad$
$\displaystyle l=1,\ldots,L$ (2)

where $\mathbf{h}_{L}$ is the output of the whole Transformer stack as well as
the input of $\mathrm{FCN}_{2}$.

The Dimension Reduction Block. Note that, the Transformer does not change the
data space and thus $\mathbf{h}_{l}\in\mathbb{R}^{256\times\tau},l=0,...,L$.
The $\mathrm{FCN}_{2}$ is designed to reduce the dimensionality of
$\mathbf{h}_{L}$ so that the LOB can be represented as the latent vector
$\mathbf{Z}\in\mathbb{R}^{1\times\tilde{\tau}}$. For that purpose, we need to
first flatten $\mathbf{h}_{L}$ as a vector using 1 linear layer. Notice that,
direct concatenating each of 256 rows of $\mathbf{h}_{L}$ as a $1\times
256\tau$ vector will lead to too large input for the successive layers of
dimension reduction. Hence, we first project
$\mathbf{h}_{L}\in\mathbb{R}^{256\times\tau}$ to
$\mathbf{h}_{L+1}\in\mathbb{R}^{40\times\tau}$ with 1 linear layer. Then the
$\mathbf{h}_{L+1}$ is flattened by concatenating its 40 rows to produce a
vector $\mathbf{h}_{L+2}\in\mathbb{R}^{1\times 40\tau}$ using 1 linear layer.
At last, 3 linear layers are carried out to reduce $\mathbf{h}_{L+2}$ to
$\mathbf{Z}=\\{z(t)\\}^{\tilde{\tau}}_{t=1}$.

The Decoder. The architecture of the decoder keeps symmetric to the encoder.
The latent vector $\mathbf{Z}$ first passes a network upside-down of
$\mathrm{FCN_{2}}$ to increase its dimension back to
$\mathbb{R}^{256\times\tau}$. Subsequently, the output of $\mathrm{FCN_{2}}$
undergoes a stack of $L$ vanilla Transformers. Finally, the output of the
Transformers stack goes through another FCN that is upside-down of
$\mathrm{FCN_{1}}$ to obtain the reconstructed data
$\mathbf{X}^{r}_{\tau}=\\{\mathbf{x}^{r}(t)\\}^{\tau}_{t=1}=g(f(\hat{\mathbf{X}}_{\tau}))\in\mathbb{R}^{10\times
4\times\tau}$.

Implementation Details. The proposed autocoder is named as SimLOB. To keep the
size of SimLOB computationally tractable, the observed $T$ steps LOB is first
split into multiple segments, each of which contains $\tau=100$ time steps as
suggested by [22]. Through sensitive analysis in Appendix E, we set $L=2$ and
$\tilde{\tau}=128$ by default. Given a pair of an original
$\hat{\mathbf{X}}_{\tau}$ and its reconstructed $\mathbf{X}^{r}_{\tau}$ with
the length of $\tau=100$ steps, the reconstruction error calculates
$\mathbf{Err}_{r}=\frac{1}{4000}\sum^{10}_{i}\sum^{4}_{j}\sum^{100}_{t}(\hat{\mathbf{x}}_{i,j}(t)-\mathbf{x}^{r}_{i,j}(t))^{2}$.
And $\mathbf{Loss}_{r}=\frac{\sum^{M}_{m=1}\mathbf{Err}^{m}_{r}}{M}$ defines
the training loss given $M$ pairs of training data sequences, by averaging the
reconstruction errors on the training data with batch size $M=128$. SimLOB is
trained by Adam [34] with learning rate 1E-4 for 200 epochs.

##  4 Empirical Studies

The following three research questions (RQs) are majorly concerned in this
work.

RQ-1: Can we learn to represent various LOB into vectors? What is the
recommended architecture?

RQ-2: Does "better representation of LOB, better calibration of FMS" generally
hold?

RQ-3: Is it really beneficial to calibrate FMS with LOB instead of traditional
mid-price?

Three groups of experiments are conducted accordingly. Group-1: the proposed
SimLOB is trained with large volume of synthetic LOB data, together with 8
state-of-the-art (SOTA) networks. The reconstruction errors on both synthetic
and real testing data are measured to assess the quality of the vectorized
representations. The parameter sensitivity of SimLOB is also analyzed in terms
of both reconstruction and calibration. Group-2: for the above 9 networks, the
represented latent vectors are applied to 10 FMS calibration tasks to show the
consistency between the reconstruction errors and the calibration errors.
Group-3: the calibration errors in terms of the vector representation learned
by SimLOB, the raw LOB and the extracted mid-price are compared.

###  4.1 General Experimental Setup

The FMS Model. The widely studied Preis-Golke-Paul-Schneid (PGPS) model is
adopted as the FMS model [35], which models all the traders as two types of
agents: 1) 125 liquidity providers who submit limited orders; 2) 125 liquidity
takers who submit market orders and cancel untraded limited orders in the LOB.
The detailed workflow and configurations can be found in Appendix B.1. In
short, the PGPS model contains 6 parameters to be calibrated.

The Dataset for Representation Learning. In the literature, augmenting the
training data with the synthetic LOB data generated by FMS is increasingly
popular [23, 27]. Hence, this work trains the SimLOB with purely the synthetic
LOB data generated by the PGPS model with different parameters. In general,
each setting of the 6-parameter tuple of PGPS actually defines a specific
market scenario. To ensure the synthetic data enjoys enough diversity, we
uniformly randomly sample 2000 different settings of those 6-parameter tuples
from a predefined range suggested by [19] (see Appendix B.2). The PGPS runs
with each setting to generate one simulated LOB for 50000 time steps,
simulating 50000 seconds of the market. Then we split each of the 2000 LOB
data as 500 sequences with 100 time steps. Thus, there are 1 million LOB
sequences with $\tau=100$, where 80% of those synthetic data is randomly
sampled for training, and the rest 20% is taken as testing data. The real
market data from the sz.000001 stock of Chinese market (during May of 2019) is
also used to test the trained networks.

The Compared Networks. In the literature, several networks have been designed
for analyzing LOB data [22]. Although their tasks main focus on a different
task of price trend prediction, they are indeed attempts of dealing with time
series LOB. This experiment selects 7 SOTA networks by taking their
architecture before the last linear layer as their representation learning
layers, and all share $\tilde{\tau}=128$. We follow the conventions of [22] to
name the 8 networks as MLP[36], LSTM[36], CNN1[37], CNN2[38], CNN-LSTM[38],
DeepLOB[32], TransLOB[30]. As their names suggested, despite that MLP, LSTM,
CNN1, and CNN2 adopt single type networks, the other works mostly employ a
Convolution-Recurrent based architecture to intuitively first extract the
prices precedence and then capture the temporal auto-correlation. TransLOB
employs the Transformer architecture to replace the recurrent network while
still keeping the convolution layers. Their input formats and sizes follow the
suggestions of their original papers and are shown in Appendix C. TransLOB-L
is an enlarged version of TransLOB to keep the size approximately aligned with
SimLOB, by using 7 stacked Transformers instead of 2 as in original TransLOB.
All the 8 compared networks are trained in the same protocol with SimLOB.

The Calibration Tasks. The 10 synthetic LOB data are generated with 10
different settings of 6-parameter tuples (see Appendix B.2), each of which
contains $T=3600$ time steps, i.e., simulating 1 hour of the market at the
frequency of 1 second. These 10 synthetic data are utilized as the target data
for calibrating PGPS. Note that these data instances are challenging that they
have a much higher frequency than traditional FMS works who can only calibrate
to daily data [19]. And the length of the target data is also at least
$10\times$ longer than the existing calibration works [10]. The objective
function for calibration employs MSE. For calibrating the mid-prices as
traditional works do, Eq.(3) gives

$\displaystyle
D(\hat{\mathbf{X}}_{T},M(\mathbf{w}))=\frac{\sum^{\lceil\frac{T}{\tau}\rceil}_{i=1}\sum^{\tau}_{t=1}(\hat{mp}((i-1)\tau+t)-mp^{\mathbf{w}}((i-1)\tau+t))^{2}}{\lceil\frac{T}{\tau}\rceil},$
(3)

where $M(\mathbf{w})=\mathbf{X}_{t}=\\{\mathbf{x}(t)\\}^{T}_{t=1}$. For
calibrating the learned latent vectors as this paper proposes, the calibration
function is slightly different. Both the target and simulated LOB with $T$
time steps is represented as $\lceil\frac{T}{\tau}\rceil$ latent vectors, each
of which has the length of $\tilde{\tau}$. Then we have Eq.(4)

$\displaystyle
D(\hat{\mathbf{X}}_{T},M(\mathbf{w}))=\frac{\sum^{\lceil\frac{T}{\tau}\rceil}_{i=1}\sum^{\tilde{\tau}}_{t=1}(\hat{z}((i-1)\tilde{\tau}+t)-z^{\mathbf{w}}((i-1)\tilde{\tau}+t)))^{2}}{\lceil\frac{T}{\tau}\rceil},$
(4)

where $\tau=100$ and $\tilde{\tau}=128$. The calibration problem can be
generally defined as $\min_{\mathbf{w}}D(\hat{\mathbf{X}}_{T},M(\mathbf{w}))$.
For simplicity, a standard Particle Swarm Optimizer (PSO) is employed to solve
it [39], and the detailed settings of PSO is given in Appendix C.1.

###  4.2 Group-1: Learning the Vectorized Representations of LOB

The reconstruction errors on 0.2 million synthetic testing data of all 9
compared algorithms are depicted in Figure 3(a). Among them, SimLOB not only
achieves the smallest averaged reconstruction error, but performs the
stablest. Specifically, the distribution of the reconstruction errors of
SimLOB follows a power law distribution with a mode of 0.003 (See Figure 7 of
Appendix D). Almost all the CNN-based networks perform the worst, suggesting
that CNN is not effective as expected to deal with LOB data. This is quite
contradictory to the trend of existing works on deep learning for LOB.

(a) The reconstruction errors on 0.2M testing LOB.

(b) The calibration errors on 10 target LOB.

Figure 3: The boxplot of the performance of 9 compared networks. The
upper/lower bound of the box indicates the top 75%/25% errors, and the red
line means the median. By definition, the point larger than 1.5$\times$ the
box length beyond the 75% bound is regarded as outlier and omitted.

Next, as shown in the first 4 rows of Figure 4, the 10 levels of bid prices of
4 randomly selected synthetic data and the reconstructed LOB of 4 networks are
visualized. For the visualizations on more data and networks, please refer to
Figure 8 of Appendix D. It is clear that the reconstructed data of SimLOB
successfully captures many details of the target LOB, especially the
fluctuations and the precedence of the 10 price levels. Comparatively, the
other networks perform quite poorer as their reconstructed price levels are
smoothed. This implies that they do not preserve the information between both
time steps and price levels.

Note that, these network are trained purely on synthetic data. We directly
apply them to reconstruct real market data. As seen in the last 4 rows of
Figure 4, the real LOB fluctuates quite slightly comparing to synthetic data.
SimLOB can still capture the information of price precedence well and the
fluctuations to some extent, though the range of prices deviates. Contrarily,
the compared networks all fail to reconstruct meaningful details. This
suggests that SimLOB can be further generalized to more diverse LOB with
special treatments, such as training with real data. The structural settings
of SimLOB are analyzed in Appendix E, which suggests $L=2$ and
$\tilde{\tau}=128$.

Figure 4: The visualization of the reconstructed data on 4 synthetic LOB data
and 4 real market LOB data. Only the 10 bid prices of each LOB is depicted
here for clarity. The $\mathbf{Err}_{r}$ between the constructed data and
target data is shown in the title of each figure.

###  4.3 Group-2: Better Representation, Better Calibration

Each of the 10 synthetic LOB data with $T=3600$ is first applied to 9 compared
networks to obtain the latent vectors, respectively. Then we apply Eq.(4) and
PSO to PGPS to obtain the best-found parameter tuple and the corresponding
simulated data. For each network, the calibration error between the target
data and the simulated data is calculated using $\mathbf{Err}_{r}$. The
detailed calibration errors are listed in Table 1. The best result between
TransLOB and TransLOB-L is listed here. Among all the learned representations,
the ones learned by SimLOB lead to the best calibration performance on 7 out
of 10 instances and the runner-up performance on data 1 and data 9. The
distribution of the calibration errors (outliers omitted) in Figure 3(b) show
the advantages of SimLOB more clearly. Furthermore, by comparing Figure 3(a)
and Figure 3(b), it can be observed that the reconstruction errors and the
calibration errors are consistent, which implies that the better the
representation of LOB is, the better the calibration of FMS is highly likely
to be.

Table 1: The calibration errors of 8 compared networks on 10 synthetic data. |  SimLOB |  MLP |  CNN1 |  CNN2 |  LSTM |  CNNLSTM |  DeepLOB |  TransLOB
---|---|---|---|---|---|---|---|---
data1 |  1.27E03 |  2.41E03 |  1.69E03 |  1.65E03 |  8.98E02 |  2.32E03 |  9.60E03 |  3.53E03
data2 |  1.30E03 |  1.85E03 |  1.61E03 |  5.67E03 |  2.43E04 |  1.51E04 |  4.14E03 |  2.52E04
data3 |  4.35E02 |  6.95E02 |  6.05E02 |  6.57E03 |  1.13E03 |  1.47E03 |  2.18E03 |  2.55E03
data4 |  5.51E02 |  1.10E03 |  2.79E03 |  1.53E03 |  4.57E03 |  8.27E02 |  2.27E03 |  1.15E04
data5 |  1.20E03 |  8.70E03 |  1.19E04 |  1.74E03 |  2.49E04 |  4.72E03 |  1.63E04 |  4.77E03
data6 |  8.58E03 |  2.11E04 |  2.87E04 |  1.34E04 |  1.07E05 |  3.40E04 |  8.42E04 |  3.29E04
data7 |  1.49E06 |  3.07E05 |  1.11E06 |  1.09E06 |  3.74E06 |  1.62E06 |  1.92E06 |  1.60E06
data8 |  1.07E03 |  2.38E03 |  3.03E03 |  1.84E03 |  1.37E04 |  2.61E03 |  8.80E03 |  5.08E03
data9 |  1.13E04 |  1.97E04 |  8.42E03 |  1.46E04 |  6.68E04 |  2.21E04 |  6.41E04 |  2.68E04
data10 |  9.83E02 |  1.36E03 |  1.98E03 |  4.30E04 |  1.31E04 |  8.28E03 |  1.51E04 |  4.12E03

###  4.4 Group-3: Calibrating Vectorized LOB is Beneficial

To demonstrate that calibrating LOB is beneficial to FMS, the above 10
synthetic LOB data is also calibrated using traditional objective function
Eq.(3) with merely mid-prices, denoted as Cali-midprice. Furthermore, one may
concern why not directly use the reconstruction error $\mathbf{Err}_{r}$ as
the calibration function to calibrate the raw LOB. Hence, we also calibrate
PGPS to the raw LOB, denoted as Cali-rawLOB.

Table 7 in Appendix F shows that SimLOB almost dominates Cali-midprice and
Cali-rawLOB. The simulated mid-prices of SimLOB, Cali-rawLOB, and Cali-
midprice on 10 target data are also depicted in Figure 5. It can be
intuitively seen that the simulated mid-price of SimLOB resembles the target
data much more than Cali-rawLOB and Cali-midprice. Comparing SimLOB with Cali-
rawLOB, $\mathbf{Err}_{r}$ essentially measures the MSE value between two very
long vectors with the length of $40\tau=4000$ on each segmentation. Though it
can tell the difference between two raw LOB data, it does not reflect the
importance of the key properties of temporal auto-correlation and the price
precedence. Comparing to Cali-midprice, SimLOB helps PGPS achieve much better
performance on not only the whole LOB but also the mid-price, simply because
the mid-price is basically a derivate of LOB. To summarize, calibrating with
more information of the market can help improve the fidelity of the
simulation, but needs effective latent representations. This supports the
initial motivation of this work.

Figure 5: The mid-prices of the target data, SimLOB, Cali-rawLOB, and Cali-
midprice.

##  5 Conclusions

This paper proposes to learn the vectorized representations of the LOB data,
the fundamental data in financial market. The autoencoder framework is
employed for this purpose, and the latent vector is taken as the
representation. In the literature, there have been increasing research works
focusing on analyzing LOB data with neural networks. However, they did not use
the autoencoder framework, neither explicitly learned the vectorized
representations. This work studies 8 SOTA neural architectures for LOB, and
discusses that the commonly used convolution layers are not effective for
preserving the specific properties of LOB. Based on that, a novel neural
architecture is proposed with three components: the first fully connected
network for extracting the features of price precedence from LOB, the stacked
Transformers for capturing the temporal auto-correlation, and the second fully
connected network for reducing the dimensionality of the latent vector.
Empirical studies verify the advantages of the proposed network in terms of
the reconstruction errors against existing neural networks for LOB. Moreover,
it is found that the effectiveness of the representation learning positively
correlated with the downstream tasks of calibrating financial market
simulation model. Further empirical findings support that the financial market
simulation models should be calibrated with compactly represented LOB data
rather than merely mid-prices or raw LOB data.

## References

  * [1] A. J. Menkveld and B. Z. Yueshen, “The flash crash: A cautionary tale about highly fragmented markets,” Management Science, vol. 65, no. 10, pp. 4470–4488, 2019.
  * [2] J. Paulin, A. Calinescu, and M. Wooldridge, “Understanding flash crash contagion and systemic risk: A micro–macro agent-based approach,” Journal of Economic Dynamics and Control, vol. 100, pp. 200–229, 2019.
  * [3] X. Wang and M. P. Wellman, “Spoofing the limit order book: An agent-based model,” in Proceedings of the 16th International Conference on Autonomous Agents and Multiagent Systems (AAMAS 2017), (Brazil), pp. 651–659, May 2017.
  * [4] J. D. Farmer and D. Foley, “The economy needs agent-based modelling,” Nature, vol. 460, no. 7256, pp. 685–686, 2009.
  * [5] V. Darley, A NASDAQ market simulation: insights on a major market from the science of complex adaptive systems, vol. 1.  World Scientific, 2007.
  * [6] A. Madhavan, “Market microstructure: A survey,” Journal of financial markets, vol. 3, no. 3, pp. 205–258, 2000.
  * [7] J. M. Hofman, D. J. Watts, S. Athey, F. Garip, T. L. Griffiths, J. Kleinberg, H. Margetts, S. Mullainathan, M. J. Salganik, S. Vazire, et al., “Integrating explanation and prediction in computational social science,” Nature, vol. 595, no. 7866, pp. 181–188, 2021.
  * [8] R. L. Axtell and J. D. Farmer, “Agent-based modeling in economics and finance: Past, present, and future,” Journal of Economic Literature, 2022.
  * [9] A. Fabretti, “On the problem of calibrating an agent based model for financial markets,” Journal of Economic Interaction and Coordination, vol. 8, pp. 277–293, 2013.
  * [10] D. Platt, “A comparison of economic agent-based model calibration methods,” Journal of Economic Dynamics and Control, vol. 113, p. 103859, 2020.
  * [11] Y. Bai, H. Lam, T. Balch, and S. Vyetrenko, “Efficient calibration of multi-agent simulation models from output series with bayesian optimization,” in Proceedings of the Third ACM International Conference on AI in Finance, pp. 437–445, 2022.
  * [12] J. Grazzini, M. G. Richiardi, and M. Tsionas, “Bayesian estimation of agent-based models,” Journal of Economic Dynamics and Control, vol. 77, pp. 26–47, 2017.
  * [13] R. Kozhan and M. Salmon, “The information content of a limit order book: The case of an fx market,” Journal of Financial Markets, vol. 15, no. 1, pp. 1–28, 2012.
  * [14] C. Chiarella and G. Iori, “A simulation analysis of the microstructure of double auction markets,” Quantitative finance, vol. 2, no. 5, p. 346, 2002.
  * [15] M. C. Recchioni, G. Tedeschi, and M. Gallegati, “A calibration procedure for analyzing stock price dynamics in an agent-based framework,” Journal of Economic Dynamics and Control, vol. 60, pp. 1–25, 2015.
  * [16] J. Grazzini and M. Richiardi, “Estimation of ergodic agent-based models by simulated minimum distance,” Journal of Economic Dynamics and Control, vol. 51, pp. 148–165, 2015.
  * [17] I. Palit, S. Phelps, and W. L. Ng, “Can a zero-intelligence plus model explain the stylized facts of financial time series data?,” in Proceedings of the 11th International Conference on Autonomous Agents and Multiagent Systems-Volume 2, pp. 653–660, Citeseer, 2012.
  * [18] S. Vyetrenko, D. Byrd, N. Petosa, M. Mahfouz, D. Dervovic, M. Veloso, and T. Balch, “Get real: Realism metrics for robust limit order book market simulations,” in Proceedings of the First ACM International Conference on AI in Finance, pp. 1–8, 2020.
  * [19] D. Platt and T. Gebbie, “Can agent-based models probe market microstructure?,” Physica A: Statistical Mechanics and its Applications, vol. 503, pp. 1092–1106, 2018.
  * [20] S. Barde, “A practical, accurate, information criterion for nth order markov processes,” Computational Economics, vol. 50, no. 2, pp. 281–324, 2017.
  * [21] F. Lamperti, “An information theoretic criterion for empirical validation of simulation models,” Econometrics and Statistics, vol. 5, pp. 83–106, 2018.
  * [22] M. Prata, G. Masi, L. Berti, V. Arrigoni, A. Coletta, I. Cannistraci, S. Vyetrenko, P. Velardi, and N. Bartolini, “Lob-based deep learning models for stock price trend prediction: a benchmark study,” Artificial Intelligence Review, vol. 57, no. 5, pp. 1–45, 2024.
  * [23] S. A. Assefa, D. Dervovic, M. Mahfouz, R. E. Tillman, P. Reddy, and M. Veloso, “Generating synthetic data in finance: opportunities, challenges and pitfalls,” in Proceedings of the First ACM International Conference on AI in Finance, pp. 1–8, 2020.
  * [24] A. Coletta, M. Prata, M. Conti, E. Mercanti, N. Bartolini, A. Moulin, S. Vyetrenko, and T. Balch, “Towards realistic market simulations: a generative adversarial networks approach,” in Proceedings of the Second ACM International Conference on AI in Finance, pp. 1–9, 2021.
  * [25] A. Coletta, A. Moulin, S. Vyetrenko, and T. Balch, “Learning to simulate realistic limit order book markets from data as a world agent,” in Proceedings of the Third ACM International Conference on AI in Finance, pp. 428–436, 2022.
  * [26] J. Li, X. Wang, Y. Lin, A. Sinha, and M. Wellman, “Generating realistic stock market order streams,” in Proceedings of the AAAI Conference on Artificial Intelligence, vol. 34, pp. 727–734, 2020.
  * [27] D. Cao, Y. El-Laham, L. Trinh, S. Vyetrenko, and Y. Liu, “A synthetic limit order book dataset for benchmarking forecasting algorithms under distributional shift,” in NeurIPS 2022 Workshop on Distribution Shifts: Connecting Methods and Applications, 2022.
  * [28] L. Hamill and N. Gilbert, Agent-based modelling in economics.  John Wiley & Sons, 2015.
  * [29] P. Yang, L. Fu, J. Zhang, and G. Li, “Ocet: One-dimensional convolution embedding transformer for stock trend prediction,” in International Conference on Bio-Inspired Computing: Theories and Applications, pp. 370–384, Springer, 2022.
  * [30] J. Wallbridge, “Transformers for limit order books,” arXiv preprint arXiv:2003.00130, 2020.
  * [31] A. Ntakaris, M. Magris, J. Kanniainen, M. Gabbouj, and A. Iosifidis, “Benchmark dataset for mid-price forecasting of limit order book data with machine learning methods,” Journal of Forecasting, vol. 37, no. 8, pp. 852–866, 2018.
  * [32] Z. Zhang, S. Zohren, and S. Roberts, “Deeplob: Deep convolutional neural networks for limit order books,” IEEE Transactions on Signal Processing, vol. 67, no. 11, pp. 3001–3012, 2019.
  * [33] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, “Attention is all you need,” Advances in neural information processing systems, vol. 30, 2017.
  * [34] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014.
  * [35] T. Preis, S. Golke, W. Paul, and J. J. Schneider, “Multi-agent-based order book model of financial markets,” Europhysics Letters, vol. 75, no. 3, p. 510, 2006.
  * [36] A. Tsantekidis, N. Passalis, A. Tefas, J. Kanniainen, M. Gabbouj, and A. Iosifidis, “Using deep learning to detect price change indications in financial markets,” in 2017 25th European signal processing conference (EUSIPCO), pp. 2511–2515, IEEE, 2017.
  * [37] A. Tsantekidis, N. Passalis, A. Tefas, J. Kanniainen, M. Gabbouj, and A. Iosifidis, “Forecasting stock prices from the limit order book using convolutional neural networks,” in 2017 IEEE 19th conference on business informatics (CBI), vol. 1, pp. 7–12, IEEE, 2017.
  * [38] A. Tsantekidis, N. Passalis, A. Tefas, J. Kanniainen, M. Gabbouj, and A. Iosifidis, “Using deep learning for price prediction by exploiting stationary limit order book features,” Applied Soft Computing, vol. 93, p. 106401, 2020.
  * [39] Y. Shi and R. Eberhart, “A modified particle swarm optimizer,” in Proceedings of 1998 IEEE international conference on evolutionary computation, pp. 69–73, IEEE, 1998.
  * [40] P. Belcak, J.-P. Calliess, and S. Zohren, “Fast agent-based simulation framework with applications to reinforcement learning and the study of trading latency effects,” in International Workshop on Multi-Agent Systems and Agent-Based Simulation, pp. 42–56, Springer, 2021.

##  Appendix A Computational Resources

The experiments run on a server with 250GB memory, the Intel(R) Xeon(R) Gold
5218R CPU @ 2.10GHz with 20 physical cores, and 2 NVIDIA RTX A6000 GPU. The
training of each network on average requires around 50 hours in the data
parallel manner with 2 A6000 GPU. The calibration of each LOB data sequence
costs 2-3 hours with 20 CPU cores or say 40 threads. The simulator can only
run on a single core, while the PSO can be run in parallel with 40 simulators.

##  Appendix B The PGPS Financial Market Simulation Model

###  B.1 The Details of PGPS

The Preis-Golke-Paul-Schneid (PGPS) model is a representative FMS model [35],
which models all the traders in the market as two types of agents: 125
liquidity providers and 125 liquidity takers. The entire simulation flowchart
can be seen in Figure 6.

Figure 6: The entire simulation workflow of the PGPS model.

The detailed settings are as follows. At each $t$-th time step, each liquidity
provider may submit a limited order with a fixed probability of $\alpha$. The
probability of each limited order being either the bid side or ask side is set
to 0.5. Each liquidity taker may submit a market order at a fixed probability
$\mu$ and cancel an untraded limited order with a probability $\delta$. The
probability of the market order being either the bid side or ask side is
$q_{taker}(t)$ and $1-q_{taker}(t)$, respectively. The $q_{taker}(t)$ is
determined by a mean-reverting random walk with a mean of 0.5, and the mean
reversion probability equals to $0.5+|q_{taker}(t)-0.5|$. The increment size
towards the mean is controlled by $\pm\Delta s$. Each liquidity taker may also
cancel an untraded limited order with a probability $\delta$. For each order,
the volume is fixed to 100 shares. The price of an ask limited order is
determined by $p^{a}_{1}(t)-\lfloor-\lambda(t)\log u)\rfloor-1$ and the price
of a bid limited order is calculated as $p^{b}_{1}(t)-\lfloor-\lambda(t)\log
u)\rfloor+1$, where
$\lambda(t)=\lambda_{0}(1+\frac{|q_{taker}(t)-0.5|}{\sqrt{\langle
q_{taker}-0.5\rangle^{2}}}C_{\lambda})$. Here $\langle
q_{taker}-0.5\rangle^{2}$ indicates a pre-computed value obtained by taking
the average of $10^{5}$ Monte Carlo iterations of $(q_{taker}-0.5)^{2}$ before
simulation. $u\sim U(0,1)$ is a uniform random number. The price of the market
order is automatically set to the best price level at the opposite side. In
summary, the PGPS model contains 250 agents and the key parameters to be
calibrated are
$\mathbf{w}=[\delta,\lambda_{0},C_{\lambda},\Delta_{s},\alpha,\mu]$.

###  B.2 Generating Synthetic LOB by Randomly Sampling 6 PGPS Parameters

The ranges of the parameters tuple $\mathbf{w}$ are listed in Table 2, as
suggested by [19].

Table 2: The ranges for randomly sampling parameters $\mathbf{w}$ to generate the synthetic data. Para. | Range | Remarks
---|---|---
$\lambda_{0}$ | $[1,200]$ | Controlling the price of each limited order.
$C_{\lambda}$ | $[1,20]$ | Controlling the price of each limited order.
$\Delta_{S}$ | $[0.0005,0.003]$ | Controlling the probability of the side of a market order (ask or bid).
$\alpha$ | $[0.05,0.45]$ | The probability for each liquidity provider to submit a limited order.
$\mu$ | $[0.005,0.085]$ | The probability for each liquidity taker to submit a market order.
$\delta$ | $[0.005,0.05]$ | The probability for each liquidity taker to cancel an untraded order.

The 10 calibration tasks with respect to synthetic data are also generated by
randomly sampling from these ranges and simulating with $M(\mathbf{w})$. Here,
we list the 10 groups of parameters in Table 3 for clarity.

Table 3: Parameter settings of generating 10 synthetic data with the PGPS model. Parameters | $\mathbf{\lambda_{0}}$ | $\mathbf{C_{\lambda}}$ | $\mathbf{\alpha}$ | $\mathbf{\mu}$ | $\mathbf{\Delta_{S}}$ | $\mathbf{\delta}$
---|---|---|---|---|---|---
data 1 | 80 | 8 | 0.1 | 0.02 | 0.002 | 0.02
data 2 | 120 | 11 | 0.2 | 0.03 | 0.003 | 0.03
data 3 | 130 | 12 | 0.3 | 0.04 | 0.003 | 0.04
data 4 | 90 | 9 | 0.15 | 0.02 | 0.001 | 0.02
data 5 | 70 | 7 | 0.15 | 0.015 | 0.0015 | 0.03
data 6 | 134.64 | 15.45 | 0.3275 | 0.07116 | 0.002 | 0.0324
data 7 | 17.7 | 11.36 | 0.2639 | 0.067 | 0.00066 | 0.03644
data 8 | 153.53 | 9.27 | 0.2983 | 0.07343 | 0.00238 | 0.01278
data 9 | 48.13 | 3.54 | 0.4374 | 0.02645 | 0.00077 | 0.01674
data 10 | 7.13 | 7.04 | 0.1106 | 0.05609 | 0.00217 | 0.01389

##  Appendix C The Characteristics of the Compared Networks

The characteristics of the 9 compared networks are listed as in Table 4. The
second and third columns give the acceptable input formats of each network,
and the last column shows the sizes of the networks. The MLP has the largest
size of 2.2E7 weights, which can be simply considered as the enlarged version
of solely the proposed $\mathrm{FCN}_{2}$. In this regard, by comparing SimLOB
with MLP, we know that the Transformer stacks are important to SimLOB. By
comparing SimLOB and TransLOB-L, we know that more Transformer does not lead
to better results (verified later in Figure 9 of Appendix E), the CNN should
be abandoned and the linear layers of $\mathrm{FCN}_{1}$ and
$\mathrm{FCN}_{2}$ are helpful.

Table 4: The characteristics of the 9 compared networks. | The length of time series $\tau$ | The size of $\mathbf{x}(t)$ |  $\\#$trainable parameters
---|---|---|---
MLP | 100 | 40 | 2.2E7
LSTM | 100 | 40 | 1.6E4
CNN1 | 100 | 40 | 3.5E4
CNN2 | 300 | 40 | 2.8E5
CNN-LSTM | 300 | 42 | 5.3E4
DeepLOB | 100 | 40 | 1.4E5
TransLOB | 100 | 40 | 1.1E5
TransLOB-L | 100 | 40 | 1.5E7
SimLOB | 100 | 40 | 1.3E7

###  C.1 The Calibration Algorithm

Note that $M(\mathbf{w})$ is basically a software simulator. In this work, it
is programmed on the Multi-Agent Exchange Environment
(MAXE)111https://github.com/maxe-team/maxe environment developed by University
of Oxford with the MIT Lisence [40]. Thus the calibration problem does not
enjoy useful mathematical properties like gradients. It is reasonable that
more advanced black-box optimization algorithms can lead to better calibration
errors. Here, for simplicity, a standard Particle Swarm Optimizer (PSO) is
employed for optimizing the above problem [39]. The hyper-parameters of PSO
follow the suggested configurations, where the population size is set to 40,
the inertia weight is set to 0.8, the cognitive and social crossover
parameters $c1=0.5,c2=0.5$. The total iteration number of PSO in each run is
fixed to 100 [10].

##  Appendix D More Results of the Reconstruction Errors

The trained SimLOB is tested on 0.2 million LOB sequences with 100 time steps.
The distribution of the reconstruction error $\mathbf{Err}_{r}$ on each
testing data is depicted in Figure 7, which follows a power law distribution
with a mode of 0.003, implying that SimLOB performs quite stable on different
LOB data.

Figure 7: The reconstruction errors on 0.2 million testing LOB sequences
follows a power law distribution.

In addition to Figure 4 in the main context, we depict the reconstructed data
of the compared networks on 5 more synthetic data to demonstrate the
superiority of SimLOB. In Figure 8, the reconstructed data of all the compared
networks except for TransLOB and TransLOB-L are visualized. The situation
remains the same that the reconstructed data of SimLOB resembles the target
data to the most. The other networks have smoothed the target data to some
extent, indicating that they are unable to fully preserve the key properties
of LOB in their latent vectors. The reason of not depicting TransLOB and
TransLOB-L is simply that their reconstructed data look like almost flat.

Figure 8: More reconstructions on the additional synthetic data are
visualized.

##  Appendix E Sensitive Analysis on the Structural Parameters of SimLOB

To analyze the architecture settings in SimLOB, we made the following
counterparts of SimLOB. The number of stacks in Transformer block is adjusted
as 2, 4, 6 and 8, respectively. The length of the latent vector $\tilde{\tau}$
is set to 64, 128, 256, and 512, respectively. The resultant variants of
SimLOB are denoted as $\mathrm{SimLOB}_{\tilde{\tau}=64}$,
$\mathrm{SimLOB}_{\tilde{\tau}=128}$, $\mathrm{SimLOB}_{\tilde{\tau}=256}$,
and $\mathrm{SimLOB}_{\tilde{\tau}=512}$, where SimLOB with
${\tilde{\tau}=128}$ is found the best choice for default settings.

###  E.1 On the number of Transformer Stack

In the main context, the number of Transformer stack is recommended as $L=2$.
Actually, with $L=2,4,6,8$, we found that the choice of $L$ did not influence
the reconstruction errors much, as shown in Figure 9. That is, the
reconstructed errors on 0.2 million testing data are around 0.025 for the 4
variants of SimLOB, and the standard deviation remains similarly, i.e., about
0.095. As adding one Transformer stack leads to around 1.5 millions more
weights, we thus use the smallest tested value of $L=2$ as the default setting
for SimLOB.

Figure 9: The number of Transformer stack in SimLOB does not influence the
reconstruction errors much.

###  E.2 On the Length of Representation Vector

The default setting of the length of the learned representation vector is
$\tilde{\tau}=128$. Such setting is based on the sensitive analysis on the
parameter $\tilde{\tau}$ with different choices like 64, 128, 256, and 512.
That is, we train SimLOB with those settings separately and tested on the 0.2
million testing data. The reconstructed errors are listed in Table 5. It is
shown that the longer the learned representation is, the better the
reconstruction errors can be obtained. This immediately suggests using larger
$\tilde{\tau}$ for representation learning.

Table 5: The reconstructed errors of different representation vectors. Networks | Reconstructed Errors
---|---
$\mathrm{SimLOB}_{\tilde{\tau}=64}$ | $0.0498\pm 0.1200$
$\mathrm{SimLOB}_{\tilde{\tau}=128}$ | $0.0250\pm 0.0997$
$\mathrm{SimLOB}_{\tilde{\tau}=256}$ | $0.0182\pm 0.1056$
$\mathrm{SimLOB}_{\tilde{\tau}=512}$ | $0.0176\pm 0.0844$

On the other hand, if we closely look at Figure 10, this advantage decreases
quickly after $\tilde{\tau}=128$. At the same time, larger lengths will make
the downstream tasks more complex and do not necessarily lead to better
calibration errors (see Table 6), where $\tilde{\tau}=128$ leads 5 best
calibration errors out of 10 instances. Thus, by balancing the reconstruction
errors and the complexity, we choose $\tilde{\tau}=128$, which has already
outperformed the SOTA, as the default setting for SimLOB.

Figure 10: The reconstructed errors of different representation vectors. Table 6: The calibration errors of different representation vectors. | $\mathrm{SimLOB}_{\tilde{\tau}=64}$ | $\mathrm{SimLOB}_{\tilde{\tau}=128}$ | $\mathrm{SimLOB}_{\tilde{\tau}=256}$ | $\mathrm{SimLOB}_{\tilde{\tau}=512}$
---|---|---|---|---
data1 | 3.53E03 | 1.27E03 | 8.68E02 | 1.03E03
data2 | 7.71E03 | 1.30E03 | 6.88E03 | 1.25E03
data3 | 5.36E02 | 4.35E02 | 4.75E02 | 4.79E02
data4 | 1.86E03 | 5.51E02 | 5.26E02 | 4.10E02
data5 | 1.62E03 | 1.20E03 | 2.49E03 | 1.38E03
data6 | 5.90E03 | 8.58E03 | 7.23E04 | 8.39E03
data7 | 2.08E05 | 1.49E06 | 1.65E05 | 1.60E05
data8 | 1.72E03 | 1.07E03 | 1.23E03 | 1.11E03
data9 | 1.51E04 | 1.13E04 | 2.64E04 | 4.01E04
data10 | 5.82E03 | 9.83E02 | 1.35E04 | 7.39E03

##  Appendix F Comparisons among SimLOB, Cali-rawLOB, and Cali-midprice on 10
Synthetic Data

To demonstrate that calibrating LOB is beneficial to FMS, the above 10
synthetic LOB data is also calibrated using traditional objective function
Eq.(3) with merely mid-prices, denoted as Cali-midprice. Furthermore, one may
concern why not directly use the reconstruction error $\mathbf{Err}_{r}$ as
the calibration function to calibrate the raw LOB. Hence, we also calibrate
PGPS to the raw LOB, denoted as Cali-rawLOB. Table 7 shows that SimLOB almost
dominates Cali-midprice and Cali-rawLOB that SimLOB wins 9 out 10 instances.
On data 2, the result of SimLOB is also very competitive to Cali-rawLOB and
far better than Cali-midprice. For Cali-rawLOB and Cali-midprice, it is hard
to tell which is significantly better. This may be the reason why traditional
FMS works did not calibrate to raw LOB data.

Table 7: The calibration errors of SimLOB, Cali-rawLOB, and Cali-midprice on 10 synthetic data. | SimLOB | Cali-rawLOB | Cali-midprice
---|---|---|---
data1 | 1.29E03 | 1.01E04 | 7.84E03
data2 | 2.70E03 | 2.26E03 | 1.29E04
data3 | 4.35E02 | 5.56E02 | 3.04E03
data4 | 5.51E02 | 6.18E02 | 3.89E03
data5 | 1.20E03 | 6.12E03 | 5.07E03
data6 | 8.58E03 | 3.29E04 | 2.71E04
data7 | 1.49E06 | 1.81E06 | 2.01E06
data8 | 1.07E03 | 2.48E03 | 2.54E03
data9 | 1.13E04 | 3.79E04 | 1.68E04
data10 | 9.83E02 | 1.49E03 | 2.88E03


# Controlling cell-matrix traction forces by extracellular geometry

Shiladitya Banerjee1 and M. Cristina Marchetti1,2 1 Department of Physics,
Syracuse University, Syracuse New York, 13244-1130, USA  2 Syracuse
Biomaterials Institute, Syracuse University, Syracuse New York, 13244-1130,
USA<EMAIL_ADDRESS>
###### Abstract

We present a minimal continuum model of strongly adhering cells as active
contractile isotropic media and use the model to study the effect of the
geometry of the adhesion patch in controlling the spatial distribution of
traction and cellular stresses. Activity is introduced as a contractile, hence
negative, spatially homogeneous contribution to the pressure. The model shows
that patterning of adhesion regions can be used to control traction stress
distribution and yields several results consistent with experimental
observations. Specifically, the cell spread area is found to increase with
substrate stiffness and an analytic expression for the dependence is obtained
for circular cells. The correlation between the magnitude of traction stresses
and cell boundary curvature is also demonstrated and analyzed.

##  1 Introduction

Living cells actively sense and respond to the physical geometry and stiffness
of their environment, which in turn affects a variety of cellular processes,
such as growth, differentiation, morphogenesis, spreading and motility [1].
Cell-matrix adhesion is mediated by integrin complexes, referred to as focal
adhesions, that bind to specific ligands on the underlying matrix. Focal
adhesions are mechanically linked to the actomyosin cytoskeleton inside the
cell, that in turn generates contractile forces on the extracellular matrix.
The interplay between substrate stiffness, intracellular contractility and the
extracellular adhesion forces controls the cell morphology and its mechanical
behavior. For instance, cells adhering to soft substrate are generally found
to spread less and have round morphology, while cells on stiff substrates have
greater spread area with more branched shapes [2].

Powerful techniques have been developed in recent years to measure the
traction forces exerted by adherent cells on compliant substrates [3].
Traction Force Microscopy is used to probe the traction stresses exerted by
cells on continuous elastic gels. The stresses are inferred from measurements
of the displacements of fiducial markers embedded in the gel before and after
cell detachment [4, 5]. In a second technique cells plated on microfabricated
pillar arrays induce bending of the elastic micropillars. The traction forces
are then obtained from by assuming a linear Hooke’s law relation between the
measured bending and the forces [6]. These experiments have demonstrated that
the mechanical response of adherent cells is controlled by a complex interplay
of substrate stiffness and geometry, myosin activity and extracellular matrix
proteins. Adhesive micro patterning has also been used as a tool for both
controlling cell shape and study the interplay between shape and cytoskeletal
organization and architecture [7]. These studies have shown that when strongly
adhesive patterns force the cell boundary to exhibit regions of high
curvature, traction stresses tend to be concentrated in these regions, while
stress fibers develop along cell boundaries linking non-adhesive zones,
confirming the crucial role of the cytoskeletal contractility and architecture
in controlling cellular stresses and morphology [8].

The role of adhesion geometry in controlling traction force distribution has
been addressed theoretically using network models and continuum mechanical
models. While models of continuum mechanical elements coupled to bio-chemical
agents have been used before to describe the traction force distribution by
adherent cells [9], continuum minimal models inspired by thermoelasticity [10]
or active gel theory [11] have recently provided new key analytical results.
Network models of the contractile cytoskeleton have also been used to describe
the relation between force distribution and shape of adherent cells [12, 13],
including networks of Hookean springs as well as cable networks that
incorporate the asymmetry of the elastic response of biopolymers such as
filamentary actin to compression version extension, with and without the
explicit inclusion of contractility. In particular, the active cable network
reproduces the arc morphology of cell boundaries pinned by strong local
adhesions that has been seen in experiments [14]. The relationship between
cell shape and adhesion geometry has also been studied by modeling cells as
contractile films bounded by the elastic cortex [15, 16, 17]. In this paper we
consider a continuum model of cells as linear, active elastic media and
demonstrate that the introduction of activity as a _spatially homogeneous_
contractile, hence negative, contribution to the pressure is sufficient to
reproduce the spatial inhomogeneous distribution of traction and cellular
stresses observed in experiments for a number of cell geometries. An
interesting extension of our work will be to introduce nonlinearity in the
continuum model to incorporate an asymmetric response to compression and
stretching. This asymmetry, arising from the nonlinear force-extension curve
of actin filaments, is known to be important in controlling the contractile
behavior of isotropic gels [18, 19] and may alter the stress distribution in
adhering cells.

In the next section we introduce our continuum model of adherent cells as
active contractile elastic media. We then use the model to study the effect of
the geometry of the adhesion region on controlling the spatial distribution of
stresses in the cell. The model can be solved analytically for a circular
cell, where we obtain an expression for the cell spread area as a function of
substrate stiffness and show that our results compares favorably to
experiments (inset of Fig. 2). The cases of elliptical, square and triangular
cells are solved numerically. We show that the geometry of the adhesive region
strongly affect the stress distribution, with traction stresses concentrated
in regions of high curvatures or at sharp corners (Fig. 1). In section 3.3 we
provide an analytical argument that quantifies the correlation between
traction stress magnitude and curvature of the cell boundary and discuss in
section 3.4 the relative roles of shear and compressional deformations in
controlling the stress distribution. We conclude with a brief discussion.

##  2 Adherent cell as a contractile gel

We consider a stationary cell adhering to an elastic substrate via stable
focal adhesion complexes. We further assume that the cell has attained its
optimum spread area on the substrate, with an average height $h$ much thinner
than its perimeter. In mechanical equilibrium, the condition of local force-
balance translates to $\partial_{\beta}\sigma_{\alpha\beta}=0$, where
${\bm{\sigma}}$ is the three-dimensional stress tensor of the cell with greek
indices taking values $x,y$ and $z$. For a thin cellular film we average the
cellular force-balance equation over the cell thickness $h$. In-plane force
balance is given by

$\partial_{j}\sigma_{ij}+\partial_{z}\sigma_{iz}=0\;,$ (1)

with $i,j$ denoting in-plane coordinates. We assume that the top surface of
the cell is stress free, $\sigma_{iz}({\bf r}_{\perp},z=h)=0$, whereas at the
cell-substrate interface $z=0$, the cell experiences lateral traction stresses
given by $\sigma_{iz}({\bf r}_{\perp},z=0)=Yu_{i}({\bf r}_{\perp},z=0)$. Here,
$Y$ denotes the substrate rigidity parameter, representing the cell-substrate
anchoring strength, and ${\bf u}({\bf r}_{\perp},z)$ is the in-plane
deformation field of the cellular medium. The thickness-averaged force balance
equation then reads [11, 10],

$h\partial_{j}\overline{\sigma}_{ij}=Yu_{i}\;,$ (2)

where $\overline{\sigma}_{ij}({\bf
r}_{\perp})=\int_{0}^{h}(dz/h)\sigma_{ij}({\bf r}_{\perp},z)$. It is
worthwhile to mention that the assumption of in-plane traction forces is a
good approximation for fully spread stationary cells making almost zero
contact angle with the substrate. During the early stages of spreading and
migration, cells can exert appreciable out-of-plane traction forces via
rotation of focal adhesions [20]. In the following we will drop the overbear
indicating the average and refer to thickness averaged quantities throughout.
The quantity $T_{i}=Yu_{i}$ is a stress in three dimensions, i.e., a force per
unit area. It describes the in-plane traction force per unit area that the
cell exerts on the substrate. The assumption of local elastic interactions
with the substrate strictly holds on elastic substrates that are much thinner
than the lateral size of the cell [21] or on micropillar substrates [10]. The
substrate rigidity parameter $Y$ depends on the stiffness of the underlying
substrate as well as on the density $\rho_{f}$ and stiffness $k_{f}$ of focal
adhesions. For an elastic substrate of shear modulus $\mu_{s}$ and thickness
$h_{s}$, $Y$ takes the simple form [21],
$Y^{-1}=\frac{1}{k_{f}\rho_{f}}+\frac{1}{\mu_{s}/h_{s}}$.

We model the cell as an isotropic and homogeneous elastic material with
additional internal active stresses due to actomyosin contractility. The
constitutive relation for the cellular stress tensor is then given by,

$\sigma_{ij}=\frac{E}{2(1+\nu)}\left(\frac{2\nu}{1-2\nu}\bm{\nabla}\cdot{\bf
u}\
\delta_{ij}+\partial_{j}u_{i}+\partial_{i}u_{j}\right)+\sigma_{a}\delta_{ij}\;,$
(3)

where $E$ and $\nu$ denote the Young modulus and Poisson ratio of the cellular
material, respectively. Actomyosin contractility is modeled as a negative
contribution to the local pressure, corresponding to $\sigma_{a}>0$. The
assumption of linear elasticity is valid on time scales shorter than
cytoskeletal turnovers, that are indeed slowed down by strong adhesion to the
substrate. Equations (2) and (3), subject to the boundary condition
$\sigma_{ij}n_{j}|_{\Omega}=0$, wtih $\Omega$ the cell boundary and ${\bf n}$
the outward unit normal on $\Omega$, completely describe the equilibrium of an
adherent cell. As a consequence of the stress free condition at the lateral
cell boundary, the net traction force transmitted by the cell to the substrate
vanishes, i.e., $\int_{A}d^{2}{\bf r}\ Yu_{i}=\oint_{\Omega}ds\
\sigma_{ij}n_{j}=0$. It is instructive to consider two limiting cases for the
anchoring strength. When the cell is rigidly anchored onto the substrate,
corresponding to $Y\rightarrow\infty$, we find ${\bf u}=0$, defining the
reference state for elastic deformations. In our model the reference cell
shape is then dictated by the geometry of the adhesion patch, which can be
controlled in experiments by micropatterning substrates by adhesion proteins.
In contrast, when $Y\rightarrow 0$, the cell does not adhere to the substrate
and the equilibrium state is uniformly contracted state, with a density
enhancement $\delta\rho=-\bm{\nabla}\cdot{\bf
u}=\sigma_{a}(1+\nu)(1-2\nu)/E(1-\nu)$. In the following, we investigate
analytically and numerically solutions of the cell elasticity equations (2)
and (3) subject to stress-free boundary conditions in various planar
geometries.

##  3 Results

###  3.1 Spatial distribution of traction stresses is sensitive to adhesion
geometry

Figure 1: Equilibrium cell shapes for various adhesion patterns : Circle (top
left), ellipse (top right), square (bottom left) and equilateral triangle
(bottom right). The color map indicates magnitude of the traction $|{\bf
T}|=Y|{\bf u}|$, and the arrows demote the direction of the traction vectors.
The reference shapes for all the four patterns have an equal area of $1000\
\mu m^{2}$. The other parameters are: $E=1~{}{\rm kPa}$, $\nu=0.4$,
$\sigma_{a}=1~{}{\rm kPa}$, $\mu_{s}=10~{}{\rm kPa}$, $h_{s}=30~{}\mu{\rm
m}m$, $h=0.2~{}\mu{\rm m}$.

The spatial distribution of traction stresses exerted by cells on substrate
and the corresponding organization of stress and deformation inside the cell
are affected by the geometry of adhesive patterns. Using micropatterning
techniques, cell shapes can be constrained to adhere to controlled geometrical
patterns [22, 23]. In our model the shape determined by the pattern in the
limit of infinite adhesion strength provides the reference shape for the cell.
Here we investigate four reference cell shapes: circle, ellipse, square and
equilateral triangle. These are chosen to have the same reference area but
different perimeters. The case of a circular cell can be treated analytically,
as described below. For the other shapes the elasticity equations (2) and (3)
are solved numerically using the MATLAB pde toolbox. We assume the
contractility $\sigma_{a}$ to be uniform and of order of the cellular Young’s
modulus. Heatmap of traction stresses are shown in Fig. 1. In all cases the
traction stresses are concentrated at the cell periphery, irrespective of the
reference shape. The magnitude of the local traction stress is, however,
higher in regions of high curvatures or at sharp corners.

For a circular cell, Eqs. (2) and (3) can be solved analytically [10, 24].
Assuming in-plane rotational symmetry, it is convenient to use polar
coordinates $r$ and $\theta$, denoting radial and angular coordinates, and
demand that no quantity depend on $\theta$. The equation for the radial
displacement $u_{r}$ about a circular reference state of radius $R_{0}$, is
then given by

$r^{2}\partial_{r}^{2}u_{r}+r\partial_{r}u_{r}-(1+r^{2}/\ell_{p}^{2})u_{r}=0\;,$
(4)

where the penetration length $\ell_{p}$ describes the localization of traction
stresses at the cell boundary. It is given by :

$\ell_{p}^{2}=\frac{Eh(1-\nu)}{Y(1+\nu)(1-2\nu)}\;,$ (5)

and is essentially controlled by the ratio of the cell stiffness $\sim E$ to
the substrate rigidity $\sim Y$. The penetration length is short on stiff
substrates and increases with decreasing substrate rigidity. The solution of
Eq. (4) with the boundary conditions $\sigma_{rr}(r=R_{0})=0$ and
$u_{r}(r=0)=0$ is given in terms of modified Bessel functions of the first
kind as,

$u_{r}(r)=-\sigma_{a}R_{0}\left[\frac{(1+\nu)(1-2\nu)}{E(1-\nu)}\right]I_{1}(r/\ell_{p})g(R_{0}/\ell_{p})\;,$
(6)

with $g(s)=\left[sI_{0}(s)-\frac{1-2\nu}{1-\nu}I_{1}(s)\right]^{-1}$. As
anticipated, the deformation $u_{r}$ vanishes for all $r$ when
$Y\rightarrow\infty$, when the adhering circular cell is maximally spread and
has its largest undeformed radius $R_{0}$.

###  3.2 Cell spread area is sensitive to substrate stiffness and
contractility

The optimal spread area of the cell is controlled by the interplay between
cell contractility, as described by the active pressure $\sigma_{a}$, and the
traction forces on the substrate. In the case of a circular cell, where the
deformation induced by adhesion is given by Eq. (6), the steady state cell
area is given by,

$A=\pi(R_{0}+u(R_{0}))^{2}\;,$ (7)

with $R_{0}$ the reference radius corresponding to the maximal spread area
$A_{\infty}=\pi R_{0}^{2}$ attained on an infinitely rigid substrate, where
$u_{r}(r)=0$. To make contact with experiments, we investigate the ratio
$A/A_{\infty}$, the relative cell spread area, as a function of substrate
stiffness and contractility.

Figure 2: Optimal shape of a triangular cell for different values of the
active pressure $\sigma_{a}$ and the substrate shear modulus $\mu_{s}$, with
$E=1~{}{\rm Pa}$. The color map represents the magnitude of the displacement
vector $|{\bf u}|$ (proportional to the traction force) about an equilateral
triangular reference shape of area $1000~{}\mu{\rm m}^{2}$. The cell spread
area increases with increasing substrate stiffness and decreases with
increasing $\sigma_{a}$. Inset (Left) : Least-square fit of the relative cell
spread area $A/A_{\infty}$ obtained from the model using Eq. (7) (solid) to
the experimental data reported in Ref. [25] (solid red circles). The fitting
parameters are $E=911~{}{\rm Pa}$ and $\sigma_{a}=1589~{}{\rm Pa}$. Inset
(Right) : Relationship between cellular Young’s modulus $E_{c}$ and
contractility $\sigma_{a}$. Here we tune $\sigma_{a}$ to desired values and
then determine the fitting parameter $E_{c}$ using data in Ref. [25]. Other
parameters : $\nu=0.4$, $h_{s}=30~{}\mu{\rm m}$, $h=0.2~{}\mu{\rm m}$.

On stiff substrates, where $R_{0}\gg\ell_{p}$, i.e., the traction stress
extends over a length much smaller than the reference cell radius,
$u_{r}(R_{0})\simeq-\sigma_{a}\ell_{p}/B$, where the compressional modulus $B$
is given by $B=E(1-\nu)/\left[(1+\nu)(1-2\nu)\right]$. The relative spread
area then takes the simple form
$A/A_{\infty}\simeq\left(1-\frac{\sigma_{a}}{R_{0}}\sqrt{h/BY}\right)^{2}$.
Letting $Y\simeq\mu_{s}/h_{s}$, we note that increasing substrate stiffness
increases relative spread area, with $A/A_{\infty}\rightarrow 1$ as
$\mu_{s}\rightarrow\infty$, in qualitative agreement with experiments [2, 26,
25]. In contrast, increasing the contractile pressure $\sigma_{a}$ reduces the
optimal cell spread area, consistent with the experimental observation that
myosin-II activity retards cell spreading [27]. To make a quantitative
comparison with experiments, we fit Eq. (7) to experimentally reported data on
the projected area of cardiac myocytes cultured on N-cadherin coated PA gels
of varying stiffnesses [25]. Here the maximal spread area $A_{\infty}$ is
taken to be equal to the cell projected area on a glass substrate (shear
modulus $\sim$ 30 GPa), which is $\simeq 690\ \mu m^{2}$. The fit, shown in
the left inset of Fig. 2, is obtained using the active contractility
$\sigma_{a}$ and the cellular Young’s modulus $E$ as the fitting parameters. A
least-square fit gives us $E=911~{}{\rm Pa}$ and $\sigma_{a}=1589~{}{\rm Pa}$.
Although the strength of contractility is likely to depend on cell type, it is
worth highlighting that the fit value for $\sigma_{a}$ is of the same order of
magnitude as previously used in Ref. [24] to fit the measured value of the
surface tension of a colony of epithelial cells. Next, we tune the
contractility $\sigma_{a}$, which can be artificially controlled through
pharmacological interventions, and determine the corresponding best fit value
of the cellular Young’s modulus $E_{c}$. Our result (Fig. 2, right inset)
indicates a linear relationship between the cellular Young’s modulus and the
contractile stress. There are indeed experimental data available [28] that
show that the cell stiffness increases linearly with contractility for
adherent cells 111We thank the anonymous referee for suggesting this fit and
pointing out to us Ref. [28].. This suggests that our model could be used to
infer contractility from measurements of cellular stiffness. Figure 2 also
demonstrates the competing roles of contractility and adhesion in controlling
optimal cell shapes for a chosen triangular reference state. On softer
substrates the triangular cell retains its topology and contracts by an amount
proportional to $\sigma_{a}$, whereas on stiffer substrates the corners tend
to form protrusions.

###  3.3 Traction forces increase with cell boundary curvature

Figure 3: (a) Force-balance on a thin slice of cellular material at the cell
boundary. (b) Force-balance at a generic sharp corner with opening angle
$\phi$. (c) Traction stress magnitude at the cell edge as a function of the
local curvature $\kappa$ for the elliptical cell of Fig. 1.

When the boundaries of the adhesion pattern exhibits non-uniform curvature,
the traction stresses are higher at regions of high curvatures. This is seen
for example in Fig. 1 for the case of an elliptical reference shape. To
justify this we propose a simple analytical argument based on local force
balance. Consider a thin slice of cellular material at the cell periphery of
width comparable to penetration length $\ell_{p}$ and arc length
$R\Delta\theta$ much less than the cell perimeter (Fig. 3(a)), with $1/R$ the
local curvature of the cell element. At the outer edge of this element, the
only force on the cell is the reaction to the traction by the cell on the
substrate traction, of areal density $-{\bf T}$, with ${\bf T}=Y{\bf u}$. This
yields an outward total force on the outer edge of the cell element of
magnitude $TR\Delta\theta\ell_{p}$, with $T>0$. At the interior edge, the
cellular element experiences a contractile force of magnitude
$\sigma_{n}(R-\ell_{p})\Delta\theta\ell_{p}$, where $\sigma_{n}$ is the normal
stress pulling the inner contour inwards and has contributions from active as
well as passive elastic stresses. The lateral stresses $\sigma_{t}$
contributes to an effective line tension $\sigma_{t}\ell_{p}R\Delta\theta$ of
the cell element. Due to the curvature of the boundary element, the line
tension generates an inward Laplace pressure of magnitude
$\sigma_{t}\ell_{p}/R$. Local balance of forces then yields,

$TR\Delta\theta\ell_{p}-\sigma_{n}(R-\ell_{p})\Delta\theta\ell_{p}=R\Delta\theta\ell_{p}\sigma_{t}\frac{\ell_{p}}{R}\;.$
(8)

The above law can be written down in a compact form as,

$T=\sigma_{n}+(\sigma_{t}-\sigma_{n})\ell_{p}\kappa\;,$ (9)

with $\kappa=1/R$, the local curvature of the boundary element. Equation (9)
then tells us that local magnitude of traction increases linearly with
increasing boundary curvature. The lateral and normal stresses $\sigma_{t}$
and $\sigma_{n}$ can be expressed in terms of the local cellular stresses in
polar coordinates as
$\sigma_{t}=\sigma_{\theta\theta}-\partial_{\theta}\sigma_{r\theta}$ and
$\sigma_{n}=\sigma_{rr}$. The linear dependence of $T$ on $\kappa$ strictly
holds in the limit $\ell_{p}\kappa\ll 1$. In addition, non-local elastic
interactions can also affect the dependence of traction magnitude on local
curvature. Figure. 3(c) shows the dependence of the magnitude of the traction
stress at the cell boundary on local curvature for an elliptical cell as shown
in Fig. 1. For low $\kappa$, the traction stress magnitude increases linearly
with $\kappa$ before reaching a plateau at higher values of $\kappa$.

When the cell boundary exhibits a sharp corner with opening angle $\phi$, as
shown in Fig. 3(b), the local force-balance is given by,

$T=\sigma_{n}+2\sigma_{t}\cos{(\phi/2)}\;,$ (10)

where $\sigma_{n}$ acts along the bisecting line of the corner. Hence smaller
the opening angle, the larger is the traction force.

###  3.4 Mechanical anisotropy induced by geometric anisotropy

Figure 4: Cell shape anisotropy correlates with internal stress anisotropy.
(a) Heatmap of internal compressive stress $\sigma$ (left) and maximum shear
stress $\sigma_{s}$ (right) corresponding to various reference shapes :
circle, ellipse, square and equilateral triangle. The reference shapes all
have an equal area of $1000\ \mu m^{2}$. (b) Average maximum shear
$\bar{\sigma}_{s}$ as a function of eccentricity $e$ for elliptical cells of
same reference area ($1000\ \mu m^{2}$). Equilibrium shapes with colorplot of
$\mu$ are given as plot markers. Parameters : $E=1$ kPa, $\nu=0.4$,
$\sigma_{a}=1$ kPa, $\mu_{s}=10$ kPa, $h_{s}=30\ \mu m$, $h=0.2\ \mu m$.

The spatial distribution of internal stresses $\sigma_{ij}$ within the cell
depends on cell shape, which is in turn controlled by the geometry of the
adhesive region. Experimentally ${\bm{\sigma}}(x,y)$ can be obtained from the
measured distribution of traction stresses ${\bf T}(x,y)$, inverting the local
force-balance condition $\partial_{j}\sigma_{ij}=T_{i}$ [29]. The elasticity
equations Eqs. (2) and (3) can be recast as a single partial differential
equation for the internal stress tensor $\sigma_{ij}$, given by

$\ell_{p}^{2}\left[\partial_{i}\partial_{k}\sigma_{kj}\right]^{S}+\delta_{ij}\sigma_{a}=\sigma_{ij}+\frac{1-2\nu}{\nu}\delta_{ij}\left(\sigma_{kk}-2\sigma_{a}\right)\;,$
(11)

where $[...]^{S}$ denotes symmetrization with respect to indices that are not
summed over, i.e.,
$\left[\partial_{i}\partial_{k}\overline{\sigma}_{kj}\right]^{S}=\frac{1}{2}\left[\partial_{i}\partial_{k}\overline{\sigma}_{kj}+\partial_{j}\partial_{k}\overline{\sigma}_{ki}\right]$.
We have investigated numerically the solution of Eq. (11) with stress free
boundary condition $\sigma_{ij}n_{j}=0$. To understand the role of shear and
compressional deformations in different geometries, it is instructive to
diagonalize the stress tensor and display the results in terms of linear
combinations of the eigenvalues $\sigma_{1}$ and $\sigma_{2}$. The sum
$\sigma=\frac{1}{2}(\sigma_{1}+\sigma_{2})$ is simply half the trace of the
stress tensor and describes compressional deformations. The difference
$\sigma_{s}=\frac{1}{2}|\sigma_{1}-\sigma_{2}|=\sqrt{[\sigma_{xx}-\sigma_{yy}]^{2}+4\sigma_{xy}^{2}}$,
is controlled by normal stress $\sigma_{xx}-\sigma_{yy}$ and shear stress
$\sigma_{xy}$. For an isotropic reference shape, such as the circle,
$\sigma_{1}=\sigma_{2}$ and $\sigma_{s}=0$, whereas for anisotropic shapes
such as the ellipse, one expects nonzero values for the local maximum shear
$\sigma_{s}$.

Fig. 4(a) shows heatmaps of the spatial distribution of $\sigma$ and
$\sigma_{s}$ for various reference shapes - circle, ellipse, square and
equilateral triangle. Irrespective of the shape of the adhesion geometry,
$\sigma$ is maximum at the cell center, indicating build-up of compressive
stresses. The compressional stress $\sigma$ always vanishes at the boundary,
and it does so more rapidly at regions of high curvature or at sharp corners.
In contrast, the shear stress $\sigma_{s}$ is identically zero for isotropic
shapes, defined as those that have a gyration tensor that is diagonal, with
equal eigenvalues. The circle, triangle and square are all in this class.
Local stress anisotropy as measured by $\sigma_{s}$ is nonzero for elliptical
shapes and shear stresses build up at the center of the ellipse. The shape
anisotropy of ellipses can be quantified by their eccentricity
$e=\sqrt{1-(b/a)^{2}}$, with $a$ and $b$ the semi-major and semi-minor axes.
Figure 4(b) shows the spatial average of $\sigma_{s}$ over the area $A$ of the
cell, defined as $\bar{\sigma}_{s}=\frac{1}{A}\int_{A}d^{2}{\bf r}\
\sigma_{s}$, as a function of the eccentricity $e$. The average shear stress
$\bar{\sigma}_{s}$ increases with $e$ with a sharp rise as $e\rightarrow 1$,
indicating a positive relationship between geometrical and mechanical
anisotropy in adherent cells. Our theoretical model thus confirms the
experimental result that cell mechanical anisotropy increases with increasing
aspect ratio, as previously reported for single endothelial cells with the
same spread area [30].

##  4 Discussion

We have used a continuum model of an adherent cell on a substrate as an active
contractile medium to study the role of adhesion geometry in controlling cell
shape, cell spreading and the spatial distribution of traction stresses. More
realistic future modeling should take into account that a cell is a highly
heterogeneous material with spatially varying stiffness [31]. It is however
intriguing to note that the simplified assumption of homogeneity and isotropy
in the underlying cytoskeletal network can reproduce several of the known
experimental results. The central input of the model is the cell contractility
or activity $\sigma_{a}$, a negative contribution to the pressure that enters
the constitutive equation for the cellular material. In general, $\sigma_{a}$
will be determined by the concentration and activity of myosin proteins cross
linking the actin cortex and controlling the formation of stress fiber. In our
model $\sigma_{a}$ is assumed to be a constant parameter, to be determined by
fitting experiments. We consider cells adhering to flat substrates that have
been patterned with adhesive patches, consisting for instance of fibronectin
coatings, of specific geometry and examine the role of the geometry of the
adhesive patch in controlling the spatial distribution of stresses in the
cellular material. The reference state for our cell is the limit of infinitely
strong adhesion, where the cell shape and lateral extent and determined
entirely by the shape and size of the adhesive patch. For finite adhesion
strength, cell elasticity and contractility yield deviations form this
reference state. We restrict ourselves to considering continuous or densely
spaced adhesion sites. For discrete or sparsely distributed adhesion sites,
non-adherent segments in the cell boundary could likely exhibit morphological
transitions induced by contractile activity and substrate stiffness [17]. In
agreement with experimental observations, we find that cells spread more on
stiff substrates and we provide an expression for the cell area versus
substrate stiffness for the case of a circular cell. We show that this
expression fit the data for spread areas of cardiac myocytes on substrates of
various sitffness values(see inset of Fig. 2). We demonstrate analytically and
numerically that strong traction stresses correlate with regions of high cell
boundary curvature, in agreement with experimental observations. Further, as
reported in experiments on single endothelial cells, our model demonstrates
that cell mechanical anisotropy is higher on elongated cells than on rounded
ones for fixed area [30].

Understanding the relation between cell morphology, the cell’s mechanical
response and cell fate is an important question in cellular biophysics. Our
simple model highlights the correlation between the geometry of adhesion sites
and cell morphology and demonstrates that traction forces by cells can be
tuned by controlling the geometry of adhesive regions. An important open
question not addressed by this simple model where the adhesive patch geometry
solely controls the cell shape is how cell morphology is determined by the
interplay of cell-substrate adhesion and dynamical reorganization of the
cytoskeletal architecture in response to the adhesion stimulus. To understand
this it will be necessary to incorporate the dynamical feedback between actin
reorganization and adhesion kinetics.

This work was supported by the National Science Foundation through award
DMR-1004789.

## References

## References

  * [1] D.E. Discher, P. Janmey, and Y. Wang.  Tissue cells feel and respond to the stiffness of their substrate.  Science, 310(5751):1139, 2005.
  * [2] T. Yeung, P.C. Georges, L.A. Flanagan, B. Marg, M. Ortiz, M. Funaki, N. Zahir, W. Ming, V. Weaver, and P.A. Janmey.  Effects of substrate stiffness on cell morphology, cytoskeletal structure, and adhesion.  Cell Motil Cytoskel, 60(1):24–34, 2005.
  * [3] A.K. Harris, P. Wild, and D. Stopak.  Silicone rubber substrata: a new wrinkle in the study of cell locomotion.  Science, 208(4440):177, 1980.
  * [4] M. Dembo and Y.L. Wang.  Stresses at the cell-to-substrate interface during locomotion of fibroblasts.  Biophysical journal, 76(4):2307–2316, 1999.
  * [5] J.P. Butler, I.M. Tolić-Nørrelykke, B. Fabry, and J.J. Fredberg.  Traction fields, moments, and strain energy that cells exert on their surroundings.  American Journal of Physiology-Cell Physiology, 282(3):C595–C605, 2002.
  * [6] J. L. Tan, J. Tien, D. M. Pirone, D. S. Gray, K. Bhadriraju, and C. S. Chen.  Cells lying on a bed of microneedles: An apporach to isolate mechanical force.  PMNAS, 100:1484–1489, 2003.
  * [7] M. Théry.  Micropatterning as a tool to decipher cell morphogenesis and functions.  Journal of Cell Science, 123(24):4201–4213, 2010.
  * [8] A.D. Rape, W. Guo, and Y. Wang.  The regulation of traction force in relation to cell shape and focal adhesions.  Biomaterials, 32(8):2043–2051, 2011.
  * [9] V.S. Deshpande, R.M. McMeeking, and A.G. Evans.  A bio-chemo-mechanical model for cell contractility.  Proceedings of the National Academy of Sciences, 103(38):14015–14020, 2006.
  * [10] C.M. Edwards and U.S. Schwarz.  Force localization in contracting cell layers.  Physical Review Letters, 107(12):128101, 2011.
  * [11] S. Banerjee and M.C. Marchetti.  Substrate rigidity deforms and polarizes active gels.  EPL (Europhysics Letters), 96:28003, 2011.
  * [12] C.A. Lemmon and L.H. Romer.  A predictive model of cell traction forces based on cell geometry.  Biophysical journal, 99(9):L78, 2010.
  * [13] P.G. Torres, IB Bischofs, and US Schwarz.  Contractile network models for adherent cells.  Physical Review E, 85(1):011913, 2012.
  * [14] C.A. Lemmon, N.J. Sniadecki, S.A. Ruiz, J.L. Tan, L.H. Romer, and C.S. Chen.  Shear force at the cell-matrix interface: enhanced analysis for microfabricated post array detectors.  Mechanics & chemistry of biosystems: MCB, 2(1):1, 2005.
  * [15] R. Bar-Ziv, T. Tlusty, E. Moses, S.A. Safran, and A. Bershadsky.  Pearling in cells: a clue to understanding cell shape.  Proc Natl Acad Sci USA, 96(18):10140–10145, 1999.
  * [16] I.B. Bischofs, S.S. Schmidt, and U.S. Schwarz.  Effect of adhesion geometry and rigidity on cellular force distributions.  Phys Rev Lett, 103(4):48101, 2009.
  * [17] S. Banerjee and L. Giomi.  Polymorphism and bistability in adherent cells.  arXiv preprint arXiv:1209.4004, 2012.
  * [18] TB Liverpool, M.C. Marchetti, J.F. Joanny, and J. Prost.  Mechanical response of active gels.  EPL (Europhysics Letters), 85(1):18007, 2009.
  * [19] S. Banerjee, T.B. Liverpool, and M.C. Marchetti.  Generic phases of cross-linked active gels: Relaxation, oscillation and contractility.  EPL (Europhysics Letters), 96(5):58004, 2011.
  * [20] W.R. Legant, C.K. Choi, J.S. Miller, L. Shao, L. Gao, E. Betzig, and C.S. Chen.  Multidimensional traction force microscopy reveals out-of-plane rotational moments about focal adhesions.  Proceedings of the National Academy of Sciences, 110(3):881–886, 2013.
  * [21] S. Banerjee and M.C. Marchetti.  Contractile stresses in cohesive cell layers on finite-thickness substrates.  Phys Rev Lett, 109:108101, Sep 2012.
  * [22] C.S. Chen, M. Mrksich, S. Huang, G.M. Whitesides, and D.E. Ingber.  Geometric control of cell life and death.  Science, 276(5317):1425–1428, 1997.
  * [23] M. Théry, A. Pépin, E. Dressaire, Y. Chen, and M. Bornens.  Cell distribution of stress fibres in response to the geometry of the adhesive environment.  Cell Motil Cytoskel, 63(6):341–355, 2006.
  * [24] A.F. Mertz, S. Banerjee, Y. Che, G.K. German, Y. Xu, C. Hyland, M.C. Marchetti, V. Horsley, and E.R. Dufresne.  Scaling of traction forces with the size of cohesive cell colonies.  Phys Rev Lett, 108(19):198101, 2012.
  * [25] A. Chopra, E. Tabdanov, H. Patel, P.A. Janmey, and J.Y. Kresh.  Cardiac myocyte remodeling mediated by n-cadherin-dependent mechanosensing.  Am J Physiol-Heart C, 300(4):H1252–H1266, 2011.
  * [26] M. Ghibaudo, A. Saez, L. Trichet, A. Xayaphoummine, J. Browaeys, P. Silberzan, A. Buguin, and B. Ladoux.  Traction forces and rigidity sensing regulate cell functions.  Soft Matter, 4(9):1836–1843, 2008.
  * [27] T. Wakatsuki, R.B. Wysolmerski, and E.L. Elson.  Mechanics of cell spreading: role of myosin ii.  J Cell Sci, 116(8):1617–1625, 2003.
  * [28] N. Wang, I.M. Tolić-Nørrelykke, J. Chen, S.M. Mijailovich, J.P. Butler, J.J. Fredberg, and D. Stamenović.  Cell prestress. i. stiffness and prestress are closely associated in adherent contractile cells.  American Journal of Physiology-Cell Physiology, 282(3):C606–C616, 2002.
  * [29] D.T. Tambe, C.C. Hardin, T.E. Angelini, K. Rajendran, C.Y. Park, X. Serra-Picamal, E.H. Zhou, M.H. Zaman, J.P. Butler, D.A. Weitz, et al.  Collective cell guidance by cooperative intercellular forces.  Nature materials, 10(6):469–475, 2011.
  * [30] P. Roca-Cusachs, J. Alcaraz, R. Sunyer, J. Samitier, R. Farré, and D. Navajas.  Micropatterning of single endothelial cell shape reveals a tight coupling between nuclear volume in g1 and proliferation.  Biophysical journal, 94(12):4984–4995, 2008.
  * [31] S.R. Heidemann and D. Wirtz.  Towards a regional approach to cell mechanics.  Trends in cell biology, 14(4):160–166, 2004.