# PRIME: A Price-Reverting Impact Model of a cryptocurrency Exchange

Christopher J. Cho
School of Electronics and Computer Science
University of Southampton
United Kingdom
<EMAIL_ADDRESS>
Timothy J. Norman
School of Electronics and Computer Science
University of Southampton
United Kingdom
<EMAIL_ADDRESS>
Manuel Nunes
Southampton Business School
University of Southampton
United Kingdom
<EMAIL_ADDRESS>

###### Abstract

In a financial exchange, market impact is a measure of the price change of an
asset following a transaction. This is an important element of market
microstructure, which determines the behaviour of the market following a
trade. In this paper, we first provide a discussion on the market impact
observed in the BTC/USD Futures market, then we present a novel multi-agent
market simulation that can follow an underlying price series, whilst
maintaining the ability to reproduce the market impact observed in the market
in an explainable manner. This simulation of the financial exchange allows the
model to interact realistically with market participants, helping its users
better estimate market slippage as well as the knock-on consequences of their
market actions. In turn, it allows various stakeholders such as industrial
practitioners, governments and regulators to test their market hypotheses,
without deploying capital or destabilising the system.

_K_ eywords Market Impact $\cdot$ Market Microstructure $\cdot$ Market
Simulation $\cdot$ Multi-Agent System $\cdot$ Cryptocurrency

##  1 Introduction

A market simulator is an environment that allows users to test their market
strategies in a synthetic setting. These simulators are useful, as it enables
the participants to gather feedback on their ideas without deploying capital
into the markets. Consequently, many researchers and institutions commit
resources to create a realistic simulation of financial exchanges. While a lot
of work remains behind closed doors due to the proprietary nature of this
research, we have seen an increased number of open-source research in this
domain over the past few years [1, 2, 3, 4].

Most of these past works look at the exchange with a focus on the market
microstructure from a top-down perspective, looking at the behaviour of the
system as a whole. However, the bottom-up analysis is yet to be studied
extensively, especially from a transaction-by-transaction basis. This approach
is key to understanding the evolution of the limit order book over time. In
particular, research into simulations to reproduce realistic market impact
remains largely untouched, despite many studies looking at analysing and
modelling observed market impact in the real world in a theoretical manner.
This poses a problem, as a realistic market simulation should interact with
market participants in the presence of market impact and slippage. We address
this gap in research by offering two key contributions. First, we identify and
analyse the components of market microstructure observed in the BTC/USD
Futures market to be used as a benchmark for our simulation. Secondly, we
present a novel simulation framework, that leverages and enhances an existing
model, which produces desirable market microstructural dynamics in an
explainable manner.

Our contributions will have various implications for stakeholders in the world
of market simulations. Industrial practitioners will be able to test the
performance of their strategies without deploying capital into the market. On
the other hand, governments and regulators will be able to run stress tests
and scenario analyses of the markets, without intervening or destabilising the
real-world exchange. Lastly, researchers within academia can benefit from a
realistic model of an exchange to test, train and benchmark their theories
against others, in a controlled and reproducible manner.

##  2 Background Concepts

To understand the simulation of financial exchange and its usage, a basic
understanding of the exchange, as well as common methods to analyse and
simulate the exchange needs to be understood.

###  2.1 Continuous Double Auction

The financial exchange is a Continuous Double Auction (CDA), where multiple
buyers and sellers simultaneously transact a homogeneous product in a
centralised venue. These participants trade with one another by using a
standardised order format that contains information on price, volume, and the
buy/sell status. Additionally, the orders are divided into two different
types: limit orders (LO) and market orders (MO). Limit orders can be viewed as
the participant’s intention to trade. These orders are placed inside the limit
order book (LOB) at the price and quantity specified by the agent and wait for
a counterparty to the trade to appear. On the other hand, market orders are
transacted immediately against the current best limit order, by becoming the
counterparty to an existing limit order.

###  2.2 Limit Order Book

The Limit Order Book is a histogram of limit order volume at each discrete
price level of a given exchange. This histogram evolves over time, as limit
orders are added and taken out of the order book via market orders and
cancellations. This means that in theory, every outstanding order in the order
book can be mapped to a participant, alongside the timestamp of the order.
However, most exchanges limit this information by anonymising the orders to
protect the identity of their clients. In addition, there are specific
terminologies used in the Limit Order Book. A “bid” is a limit order looking
for a counterparty to buy the asset from, whereas an “ask” is a limit order
looking to sell the asset. At any point in time, the “spread” of the market is
defined as the gap between the best bid and ask of the market, and the “mid-
price” is defined as the halfway point between the best bid and ask price of
the market.

###  2.3 Market Simulator

In order to create a realistic market simulation, a simulator that can
accommodate our model is necessary. Out of a handful of open-source market
simulators, we take the ABIDES [1] (the Agent-Based Interactive Discrete Event
Simulation environment) platform and add additional utilities to accommodate
our simulation. This platform was chosen over others for a few reasons.
Firstly, it has a built-in script that allows for easy parallelisation of the
simulation workflow, allowing us to carry out large-scale empirical studies
with less effort. Secondly, the environment has built-in examples of the
“agent” class, that include simple, but useful and commonly used agents such
as Zero Intelligence and Noise Trader. Lastly, the platform allows
participating agents to query an “oracle”, that can provide the querying
agents with information that is exogenous to the simulation (e.g. true price
of the asset using a historic dataset). This last point is crucial, as it
opens up the possibilities for the simulation to track the historic price
series of an asset, which allows for a better comparison between historic and
simulated market data.

##  3 Data Collection

Financial assets are traded either on an exchange or off the exchange in an
“over the counter” (OTC) fashion. This research focuses on the former, as OTC
transactions are done privately with little information shared with the
public. Data from exchanges are in general homogeneous, well-maintained, and
plentiful. However, exchanges of traditional assets such as equities and
futures charge users to access their order book data, which is a huge barrier
to entry for academic researchers. On the other hand, data from the less
established cryptocurrency exchange comes free of charge. Therefore, our
research is carried out on the Level 3 LOB data collected from the Binance
API. More specifically, we queried data from the USD/BTC perpetual futures
market from July 2020 to October 2021. This is an order-by-order breakdown of
all executed orders during this period. Unfortunately, data is unavailable for
limit order submissions and cancellations, somewhat limiting the scope of our
research. In addition, we recognise that due to the decentralised nature of
cryptocurrencies, the data from Binance may not be exactly in-line with the
entire cryptocurrency market. However, due to their size (largest exchange by
transaction volume, March 2023) and the existence of high-frequency traders
who arbitrage between various cryptocurrency exchanges, we assume that the
data we collect from Binance is representative of the entire market for the
purpose of this research.

##  4 Market Impact in the Real World

###  4.1 What is Market Impact?

Market impact is defined as the change in the price of an asset caused by the
trading of that asset. From the perspective of a single order, the market
impact is deterministic based on the size and type of the order, and the
prevailing liquidity in the limit order book at the point in time of the
order. This means the price move is larger if the traded volume is large or
prevailing liquidity is thin. In addition, the available liquidity at each
price increases roughly linearly as the price moves further away from the best
bid/ask price (until a certain distance). This leads to the phenomenon where
price requires quadratically increasing trade volume to move linearly away
from the mid-price. This is one of the fundamental rational behind the Square
Root law [5] of market impact.

The market impact becomes less mechanical when approached from a non-
instantaneous perspective. This is due to the concept of impact reversion and
the hidden liquidity arising from what is commonly known as the Latent Order
Book [6]. This is the idea that most market participants do not place limit
orders at their private valuation until their valuation becomes (or close to)
the most competitive price, in order to increase their likelihood of a
transaction. This allows participants to withhold private information from the
market, allowing them to better utilise their information advantage. For
example, an individual willing to sell BTC at USD 100 may not place their
limit order when the best ask price in the market is at USD 98. However, when
a large buy order eats up all the liquidity between USD 98 and USD 100,
pushing up the best ask price to USD 101, this individual will place a limit
order asking for USD 100. This behaviour is encouraged, as placing a large
limit order at USD 100 when the market is at USD 98 would reveal to the market
that a large seller is in play, thus influencing participants’ decisions (e.g.
would-be buyers may be more inclined to wait until the price falls, rather
than transacting immediately). Consequently, it is only when the price level
changes, the hidden liquidity is brought out into the market, and causes a
reversion in market impact over time. The Propagator model [7] is a widely
used model that explains the observed part of this phenomenon by dividing up
the market impact into initial impacts and the reversion of these market
impacts when estimating the total market impact over the current period. We
too will use this framework presented by Gatheral et al for the purpose of our
research. Note that modelling the idea of visible and latent liquidity in a
fully consistent model remains an active area of research [8], and is beyond
the scope of our project. Instead, our work focuses on the proof of concept of
a model that can replicate the general behaviour of the simulated markets to
be in line with the real world in a qualitative manner, from the lens of the
Propagator model. Thus we leave the exact empirical tuning and mathematical
derivations of the model to future work.

###  4.2 Initial Market Impact

As discussed above in Section 4.1, the market impact should be approached from
a temporal perspective to avoid drawing conclusions based purely on the
mechanical nature of the LOB. Therefore we resample the trade data into
5-second time windows, where the traded volume is defined as the net buy/sell
volume across the interval (e.g. 200 buy orders and 500 sell orders within the
interval will register as -300 orders), and price change is the difference
between the close and open mid price of the 5-second interval.

Looking at the resampled data from the BTC/USD Futures market (Figure 1, left
image), a positive relationship between trade size and price change can be
observed. However, despite the generally positive relationship, some large
discrepancies in impact between orders of similar sizes are observed. This is
likely due to the differences in prevailing market liquidity and volatility
mentioned previously. A common model to control for prevailing liquidity,
current volatility, as well as the nature of the non-instantaneous market
impact, is the aforementioned Square Root model [5]. The general form of this
model is shown in Equation 1, where the current period market impact, $I$, is
modelled as a function of the volume of the orders in the current resample
period (5 seconds), $Q$, and the ”recent“ time horizon (1 hour), $T$. For
control factors, $\sigma_{T}$ is the price volatility over time $T$, $V_{T}$
is the average traded volume over time $T$, $\delta$ the exponent of the
diminishing market impact, usually assumed to be around 0.5 (hence the Square
Root law), and a constant, $k$.

$I(Q,T)=k\sigma_{T}(\frac{Q}{V_{T}})^{\delta}$ (1)

We take an incremental approach to analyse the initial market impact. We start
by ignoring the exponent, $\delta$, and adjust the market data by prevailing
price volatility, $\sigma_{T}$, and the average volume, $V_{T}$. This yields a
plot that appears better correlated (Figure 1, right image). Note that we
calculate $\sigma_{T}$ using the standard deviation of period-by-period price
change over time $T$, and $V_{T}$ using a linearly increasing weighted average
of orders over time $T$, where volume from the most recent period is assigned
the highest weighting factor.

Figure 1: BTC/USD Futures Market Impact (5 sec resample)

Despite the pleasing aesthetics of the plot, this visualisation is slightly
misleading, due to the high density of data points (over 800,000). Therefore,
when the orders are gathered into buckets of order size and their mean market
impacts plotted on top of the scatter plot (cyan lines, Figure 2, left image),
we find the relationship between order size and price change to be non-linear.
Isolating and looking closer into the bucketed means (Figure 2, right image),
increasing resistance of price change to traded volume can be observed
clearly, suggesting that the exponent $\delta$ in Equation 1 must have a value
less than 1.

Figure 2: BTC/USD Futures Market Impact with percentile buckets (5 sec
resample)

Using the data, we numerically find the value of $\delta$ using an
optimisation routine. The results of the optimisation yield an outcome
$\delta=0.59$. When this exponent is applied to the data following Equation 1,
the data and the bucketed mean impact are transformed to that shown in Figure
3. Although the scatter plot exhibits a slight kink around zero net volume, it
nonetheless displays the relationship between trade size and price change
explained by Equation 1 linearly. We believe this empirical evidence is
sufficient in showing the validity of the Square Root law in the BTC/USD
Futures market, and therefore adequate to be used as a baseline result to
compare our simulated results.

Figure 3: BTC/USD FuturesMarket Impact with percentile buckets (5 sec
resample)

Before concluding our evidence for initial market impact, we quickly note the
slight non-linearity in the adjusted bucket plots in Figure 3. The
explainability near zero net trade size is non-linear for two reasons.
Firstly, small net orders tend to be traded in a differing liquidity
environment on either side of the book. In other words, 100 units of buy
orders followed by 100 units of sell orders will not necessarily bring the
price back to the starting point, depending on the prevailing liquidity on the
buy and sell side of the book. Second, the reversion of market impact from
previous periods will account for a larger portion of the overall market
impact near the zero mark (since the current market impact is small),
distorting the relationship between current order size and market impact.

###  4.3 Reversion in Market Impact

We now present evidence for reversion in market impact in two ways. The first
method is by showing the intertemporal conditionality between the previous
trade size and the current market impact. If there is no reversion in market
impact, there should be no intertemporal conditionality between the previous
trade size and the current market impact. Therefore, under this assumption,
the sign of the previous period volume (net buy or net sell) should not
influence the current market impact in any way. However, when the current
market impact is divided into two groups based on the buy/sell status of the
previous period, we observe a gap between the groups (Figure 4). This analysis
indicates that periods following a previous period of net buy orders on
average cause a lower market impact than those following a period of net sell
orders. This means the previous order exerts market impact pressures in the
opposite direction on the current period, evidencing reversion in impact.

Figure 4: BTC/USD Futures market impact coefficient over time (5 sec resample,
volume & volatility adjusted, including exponent)

The second approach looks beyond the first-term lag. As described by the
propagator model [7], we expect the price impact of trade to revert in a
decaying fashion, as the time between the trade and price impact observation
increases. This follows Equation 2, where the total market impact, $M(t)$,
over period $t$ is the accumulation of all initial impacts made by trades
across the period, $f(\dot{x_{s}})$, after the decay of these impacts have
been taken into consideration following the decay kernel $G(t-s)$, plus a
noise term $\epsilon$. Here, the $\dot{x_{s}}$ represents the rate of trading
across our period $t$ and is assumed to be constant (thus we assume trades are
made at a constant rate, and the integral of $\dot{x_{s}}$ over the period $t$
is equal to the sum of all trades across this period). We have shown in the
previous section that the initial impact follows the functional form
represented by the Square Root model. Therefore function $f$ can be
represented by function $I$ from Equation 1, where the recent time horizon,
$T$, is fixed to be 1 hour. This allows us to find out the functional form of
the decay function, $G(t-s)$, in the BTC/USD Futures market in an empirical
manner, by investing the market impact $M(t)$ for varying levels of $t$. Note
that we ignore the noise term $\epsilon$ and assume it to be zero in our work.

$M(t)=\int_{0}^{t}f(\dot{x_{s}})G(t-s)\,ds+\int_{0}^{t}\epsilon\,dZ_{s}$ (2)

By regressing the adjusted order size from time $t-100$ to $t$ against
adjusted price change at time $t$ for all available $t$, the resulting
coefficients show the average market impact of a unit trade on the present and
future periods (Figure 5). The coefficients show a large initial impact at
$t=0$, followed by the reversion in impact, as exhibited by small negative
coefficients at $t>0$. When the cumulative market impact of a unit order at
time $t$ is plotted over 100 periods as before (Figure 5), a decay in market
impact following a power law can be observed.

Figure 5: BTC/USD Futures market impact decay over time (5 sec resample,
volume & volatility adjusted, including exponent)

###  4.4 Meta Orders

The final section of our analysis of the BTC/USD Futures market is regarding
meta orders. In general, large orders are known to be broken down into many
smaller orders in an attempt to hide the total transaction volume and reduce
the market slippage of the transaction. This is done in various methods,
including the well-known VWAP (Volume Weighted Average Price) and TWAP (Time
Weighted Average Price) execution strategies. Breaking down a large order into
smaller chunks leads to the idea of meta orders, which spawn consecutive
orders in the same direction. This idea can be verified by checking the
autocorrelation of order signs between consecutive orders.

As shown in Figure 6, the BTC/USD Futures market exhibit a clean power-law
decay in this autocorrelation of order signs, supporting the existence of
these meta orders. This autocorrelation is also something that we look to
include in our market simulation, in an attempt to produce the desired market
reaction to a transaction.

Figure 6: BTC/USD Autocorrelation of Order Sign (Actual Data, March 2021)

##  5 Simulated market impact

In section 4, we identified and analysed the market impact observed in the
BTC/USD Futures exchange using the Square Root law and the Propagator model.
This allowed us to divide up impact into two components: the initial impact
and the decay kernel. Our objective now is to ensure that our market
simulation replicates these observations as closely as possible. For this, we
analyse an existing approach to modelling the limit order book and then create
a novel approach to improve upon the shortcomings of the existing model.

###  5.1 Santa-Fe Model

Our starting point is the Santa Fe model of the market [9], originally
developed by researchers at the Santa Fe institute. This model is created
using stochastic processes to populate the limit order book following a zero-
intelligence approach. Within the model, limit orders of fixed size arrive
following a Poisson process across a predetermined price range. The arrival
rate at each price is independent of one another and independent of the number
of outstanding orders at each level. At the same time, each of these limit
orders is cancelled following a separate Poisson process, which is again
independent of the state. Simultaneously, market orders arrive following yet
another Poisson process that takes out limit orders. In doing so, the LOB is
populated as a function of limit order arrival & cancellation rate, and market
order arrival rate.

We replicate this model in the ABIDES simulator by representing these
processes using augmented Zero-Intelligence agents. This is done by dividing
the Zero-Intelligence agents into limit order agents (Algorithm 1) and market
order agents (Algorithm 2). Although it is possible within the ABIDES platform
to create a single type of Zero-Intelligence agent that is allocated a
probability of behaving as a limit or a market order agent, this increases the
run-time of the algorithm, without making any difference in the system
behaviour. It may also cause the agents to sometimes transact against their
own limit orders, which makes little sense from an explainability perspective.

$p_{cancel}:=P$(Cancelling an existing limit order)

Generate a random number $n\sim\mathcal{U}(0,1)$

if  _$n <p_{cancel}$_ then

Cancel oldest limit order if one exists

else

Obtain private valuation, $V\sim\mathcal{U}\\{1,100\\}$

Query best market ask, $A^{\prime}$

Query best market bid, $B^{\prime}$

Calculate mid price $M^{\prime}=(A^{\prime}+B^{\prime})/2$

if  _$V <M^{\prime}$_ then

Place limit order to buy at $V$

else

Place limit order to sell at $V$

end if

end if

Algorithm 1 Modified Zero-Intelligence algorithm for Santa-Fe model (Limit
order agent)

Generate a random number $m\sim\mathcal{U}(0,1)$

if  _$m\leq 1/2$_ then

Place market order to buy

else

Place market order to sell

end if

Algorithm 2 Modified Zero-Intelligence algorithm for Santa-Fe model (Market
order agent)

Reproducing the Santa Fe model in the ABIDES simulator yields very good
results on the market impact front. The diminishing initial market impact can
be observed in Figure 8 and Figure 8, and the reversion in market impact can
also be observed in Figures 9, 10 and 11. However, it was unable to show
desirable autocorrelation in order sign to produce the idea of meta orders
(Figure 12). To address this discrepancy, we use a method created by Taranto
et al [10] and incorporate a stochastic process called DAR(p), or the Discrete
AutoRegressive process of order $p$, to the market order agent (Algorithm 3).
This additional augmentation maintains the behaviour of the initial market
impact and its decay observed in the unaugmented Santa-Fe simulation, whilst
now exhibiting order sign correlation that decays in a desirable fashion
(Figure 13). Note that we tune the parameters of the DAR(p) process in
Algorithm 3 by using a Monte-Carlo method to minimise the difference between
the parameters of the power law fit observed in the real world (Figure 6) and
the one observed in the simulation (Figure 13).

Figure 7: Santa-Fe Market Impact with percentile bucket average impact (5 sec
resample)

Figure 8: Santa-Fe Market Impact with percentile bucket average impact (1 sec
resample)

Figure 9: Santa-Fe Market Impact with percentile bucket average impact by
previous order direction (5 sec resample) Figure 10: Santa-Fe market impact
coefficient over time (5 sec resample) Figure 11: Santa-Fe market impact decay
over time (5 sec resample) Figure 12: Autocorrelation of Order Sign (Santa-Fe
Simulation) Figure 13: Santa-Fe autocorrelation of Order Sign (Simulated Data,
March 2021)

The results obtained from the Santa Fe model are very promising. However, it
doesn’t quite provide us with the functionality fully equipped for real-world
simulation. Firstly, it lacks model explainability. The orders submitted in a
stochastic process by Zero-Intelligence agents do not represent a type of
market participant. Instead, it assumes that aggregate participation in the
market follows a Poisson process. This hampers our ability to attribute market
behaviour to a specific type of participant or test the impact of adjusting
agent-specific behaviour on the overall market. Secondly, the model relies on
this aggregated behaviour to remain in equilibrium to exhibit desirable market
impact. This means that adding additional participants to the simulation may
destabilise this equilibrium and push the simulation away from desirable
behaviour. Both these shortfalls remove some of the key advantages of using an
agent-based model. Another limitation of the Santa Fe model is its inability
to follow an underlying price series. As both limit and market orders arrive
independently around the mid-price, the equilibrium either mean-reverts or
trends away from the initial equilibrium on its own accord, based on the
system’s lag-dependent volatility [11]. This is why we opt to develop a new
model that addresses these shortcomings.

Generate a series of n random numbers $m_{i}\sim\mathcal{U}(0,1)*$ to form
array $M=(m_{1},m_{2},m_{3},...,m_{n})$

Generate a “parent” index $\ell\leq n$ from $\lambda_{\ell}$

where $\lambda_{\ell}\sim\ell^{(\gamma-3)/2}$

if  _$p\leq r\sim\mathcal{U}_{[0,1]}$_ then

$m_{t}=m_{\ell}$

else

$m_{t}=1-m_{\ell}$

end if

if  _$m_{t}=1$_ then

place buy market order

else

place sell market order

end if

shift index of M by 1

$(m_{1},m_{2},m_{3},...,m_{n})\rightarrow(m_{2},m_{3},m_{4},...,m_{n+1})$

remove the last element $m_{n+1}$

add $m_{t}$ as the first element $m_{1}$

Algorithm 3 Modified Zero-Intelligence algorithm for Santa-Fe model (Market
order agent, DAR(p) process)

###  5.2 PRIME

To address the issues with Santa Fe model, we propose a novel simulation
methodology we call PRIME, or the Price-Reverting Impact Model of a
cryptocurrency Exchange. This model is able to follow an underlying price
series and replicate the desired market microstructural properties of the
cryptocurrency exchange, whilst maintaining the interactivity and the
explainability of the simulation. We achieve this by augmenting and adding to
the existing Santa Fe model in various ways.

Firstly, the limit order agents no longer “receive” a private valuation from a
predetermined price range. Instead, they calculate their private valuation by
querying the prevailing mid-price of the simulation and adding a small
perturbation drawn from a uniform distribution. Secondly, at the beginning of
the simulation, the LOB is populated in a linearly increasing manner around a
specified starting price. This allows the limit order agents to query the mid-
price and populate the limit order book from the start. Lastly, market order
agents have access to the true underlying price series via the oracle, again
with uniformly distributed perturbation. These changes make the limit-order
agents function similarly to background traders, and market-order agents
become the true fundamental agents who exclusively execute their strategies
aggressively. Lastly, we include technical agents that either trend-follow or
mean-revert in the market by submitting market orders into the system.
Following these changes, the agent configuration we use in the PRIME model is
1000 Zero-intelligence limit order agents, 30 Zero-intelligence market order
agents, 10 trend-following agents and 10 mean-reversion agents.

The PRIME model behaves very similarly to the Santa Fe model if the underlying
value of the asset is fixed. In this case, limit orders arrive independently
following a Poisson process on both sides of the prevailing mid-price, with an
independent cancellation rate attributed to each order. Market orders also
arrive following a Poisson process, and if the prevailing market equilibrium
is equal to the underlying price series, the model remains identical to the
Santa-Fe model. However, if the model mid-price is not equal to the true
price, the fundamental valuation received by market-order agents will be
skewed away from the simulation’s mid-price. This results in the market orders
arriving in an asymmetric manner (different buy/sell order arrival rate) to
push the model equilibrium back towards the true price. This forces PRIME to
be always mean-reverting, which is not the case with the Santa-Fe model.
Conversely, if the true price of the underlying asset is not fixed at the
initial price, and instead fluctuates following a predetermined path (e.g.
historical price series), then PRIME offers a simulation that follows the
underlying price series (Figure 14), whilst largely maintaining the market
mechanics of the Santa-Fe model. In addition, thanks to the market-order
agents in PRIME exerting pressure on the system to revert to the true price,
the simulation can effectively interact with other agents, as perturbations
made to the system will be equilibrated based on the true market price
observed by fundamental agents. This makes PRIME better equipped to explain
and represent real-world markets.

Figure 14: PRIME Price Series with various observation error (5 sec resample)

We assess PRIME using the same methodology as before. Firstly, we seek
evidence of the Square Root law by looking at the adjusted impact plot (Figure
15). As before, we initially apply only volume and volatility adjustment
without the exponent, $\delta$ (See Equation 1). Although there is a disjoint
in adjusted returns, due to the limited resolution of the price levels in the
PRIME model (as per Santa-Fe, limit orders are only submitted within $\pm 50$
ticks of the mid-price), the resulting output clearly shows a diminishing
impact as trade size increases. Secondly, when an exponent of value 0.74 is
applied to the data, the relationship becomes much more linear, as shown by
the cyan bucketed mean impact plots in Figure 16.

Figure 15: PRIME Market Impact with percentile buckets (5 sec resample) Figure
16: PRIME Market Impact with percentile buckets (5 sec resample)

Next, we look for evidence of reversion by dividing up orders by previous
period order sign (Figure 17, right image). The market impact following a buy
order clearly shows a lower market impact than those following a sell order,
in line with the behaviour observed in the real world (Figure 4). In addition,
evidence of reversion is also found when looking at the intertemporal impact
coefficients and their cumulative decay in Figure 18. While there is a
reversion in impact that decays over time, the shape of the decay is very
linear after the first few terms, as opposed to the diminishing decay observed
in the real world (Figure 5).

Lastly, we see a promising result regarding the autocorrelation of order sign
in the left image of Figure 19. Like the real-world data and the augmented
Santa-Fe model, we observe a decay that roughly follows a power law. This is
especially encouraging, as the PRIME model does not embed artificial
stochastic processes into the participating agents to generate this result. In
fact, the autocorrelation occurs from the technical agents within the PRIME
model. The momentum agents who place orders in the same direction as previous
orders increase the order sign autocorrelation (Figure 19, middle image),
whereas the mean reversion agents do the opposite and reduce the
autocorrelation (Figure 19, right image).

Figure 17: PRIME market impact coefficient over time (5 sec resample, volume &
volatility adjusted, including exponent) Figure 18: PRIME market impact decay
over time (5 sec resample, volume & volatility adjusted, including exponent)
Figure 19: Santa-Fe autocorrelation of Order Sign (Simulated Data, March 2021)

###  5.3 Model Evaluation

Our analysis shows that both the Santa-Fe model and the PRIME model are able
to replicate the three market microstructural properties we analysed in the
BTC/USD Futures market: initial impact, reversion of impact, and the
autocorrelation of order signs. Whilst the Santa-Fe model exhibited cleaner
and more desirable results regarding the reversion of impact and the
autocorrelation of order signs, the model lacks fundamental bottom-up
explainability, especially after adding in the DAR(p) process to the market
order agent. In addition, the Santa Fe model lacks the ability to follow an
underlying price series, which is critical in replicating and comparing real-
world markets to the simulation. On the other hand, the PRIME model offers a
novel framework that addresses these weaknesses, at the cost of slight
deviation from the desired shape of impact reversion and the decay in order
sign autocorrelation. The increased explainability of the PRIME model as well
as its ability to trace an underlying price series makes the model
exceptionally useful as a benchmark agent-based simulation of the market, that
can be used to test and train market strategies and offers researchers with a
first bottom-up agent-based simulation that replicates the market
microstructure of the BTC/USD Futures market.

##  6 Related Work

This work is an extension of many previous efforts to replicate a live market
environment. McGroaty et al [12] provide a good summary of past studies and
divide them into four categories: theoretical equilibrium models; statistical
order book models; stochastic models; and agent-based models. In addition to
these models, there has been recent progress to represent the market and the
Limit Order Book dynamics by leveraging Machine-Learning techniques.

Our attempt approaches this problem using the agent-based approach to
incorporate as much explainability as possible using a bottom-up methodology.
The ability to include and customise market participants in an intuitive
manner is a key advantage of the agent-based model, and it has been increasing
in popularity [2, 3, 4, 13] as computational resources became more abundant
over the past two decades.

While the other methods of modelling the LOB were considered alongside the
Agent-Based model, we found them to be unsuited for our work. The Equilibrium
approach yields results that are too fragile in equilibrium due to the number
of assumptions required to allow the model to equilibrate. The statistical and
machine-learning approaches provide promising results [14, 15, 16], but they
lack explainability, due to their focus on the model fit, making intuitive
adjustments difficult to add in the market. Similarly, while stochastic models
[17, 9] offer more than the statistical models in terms of explainability,
they still have the tendency to group market participants into large
categories, and thus offer less explainability when compared to the agent-
based method.

In addition, our analysis of the behaviour of the BTC/USD Futures market is an
update on the previous work done by Donier et al [18]. Their work analysed the
impact on the BTC/USD market, with a focus on meta orders, before the rise in
popularity and size of the crypto market. Our findings provide an update on
the behaviour of the cryptocurrency market post the 2016 rally using
publically available data, with an additional focus on the shape of the impact
decay.

##  7 Conclusion & Future Work

In this paper, we provide two contributions to the field of cryptocurrency and
market simulations. We first identified and evidenced market microstructure in
the BTC/USD Futures market using publically available data, by leveraging
previous theories such as the Square root law and the Propagator model.
Secondly, and more importantly, we developed a novel simulation framework that
addresses the shortcoming of existing models whilst maintaining model
explainability, by replicating the key microstructural properties of the
crypto market using heterogeneous agents. This simulation framework can be
used by researchers in various ways, from investors testing their trading
strategies, regulators looking at market fragility risks, and academics
benchmarking and comparing novel models which incorporate state-of-the-art
algorithms.

In our future work, we will look to make improvements in a few areas. First,
we would like to improve the shape of the impact reversion to match that seen
in real-world markets. Secondly, we would like to empirically tune the PRIME
simulation against observed microstructure in the BTC/USD Futures market, by
varying the ratio and the aggressiveness of participating agents. Thirdly, we
would like to incorporate the idea of latent liquidity, by augmenting traders
to incorporate a more complex execution strategy described by algorithms such
as Zero-Intelligence Plus [19]. Lastly, we would like to run an extensive
study on the stability of our simulation, by subjecting the simulation to a
wide range of external agents that interact with the environment using diverse
strategies. Doing this will improve the accuracy of PRIME, and will allow
researchers to utilise the simulation with an additional layer of confidence.

## Acknowledgments

This work was supported by EPSRC Centre for Doctoral Training in Next
Generation Computational Modelling (EP/L015382/1). We also acknowledge the use
of the IRIDIS High Performance Computing Facility, and associated support
services at the University of Southampton, in the completion of this work.

## References

  * [1] David Byrd, Maria Hybinette, and Tucker Hybinette Balch.  ABIDES: towards high-fidelity market simulation for AI research.  CoRR, abs/1904.12066, 2019.
  * [2] Svitlana Vyetrenko, David Byrd, Nick Petosa, Mahmoud Mahfouz, Danial Dervovic, Manuela Veloso, and Tucker Hybinette Balch.  Get real: Realism metrics for robust limit order book market simulations.  arXiv preprint arXiv:1912.04941, 2019.
  * [3] Christopher J Cho and Timothy J Norman.  Bit by bit: how to realistically simulate a crypto-exchange.  In Proceedings of the Second ACM International Conference on AI in Finance, pages 1–9, 2021.
  * [4] Andrea Coletta, Aymeric Moulin, Svitlana Vyetrenko, and Tucker 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, pages 428–436, 2022.
  * [5] Nicolo Torre.  Barra market impact model handbook.  BARRA Inc., Berkeley, 208, 1997.
  * [6] Sriketan Mahanti, Amrut Nashikkar, Marti Subrahmanyam, George Chacko, and Gaurav Mallik.  Latent liquidity: A new measure of liquidity, with an application to corporate bonds.  Journal of Financial Economics, 88(2):272–298, 2008.
  * [7] Jim Gatheral.  No-dynamic-arbitrage and market impact.  Quantitative finance, 10(7):749–759, 2010.
  * [8] Jonathan Donier, Julius Bonart, Iacopo Mastromatteo, and Jean-Philippe Bouchaud.  A fully consistent, minimal model for non-linear market impact, 2015.
  * [9] J Doyne Farmer, Paolo Patelli, and Ilija I Zovko.  The predictive power of zero intelligence in financial markets.  Proceedings of the National Academy of Sciences, 102(6):2254–2259, 2005.
  * [10] Damian Eduardo Taranto, Giacomo Bormetti, Jean-Philippe Bouchaud, Fabrizio Lillo, and Bence Toth.  Linear models for the impact of order flow on prices i. propagators: Transient vs. history dependent impact.  arXiv preprint arXiv:1602.02735, 2016.
  * [11] Jean-Philippe Bouchaud.  Chapter 7 - agent-based models for market impact and volatility.  In Cars Hommes and Blake LeBaron, editors, Handbook of Computational Economics, volume 4 of Handbook of Computational Economics, pages 393–436. Elsevier, 2018.
  * [12] Frank McGroarty, Ash Booth, Enrico Gerding, and VL Raju Chinthalapati.  High frequency trading strategies, market fragility and price spikes: an agent based model perspective.  Annals of Operations Research, 282(1-2):217–244, 2019.
  * [13] Xintong Wang and Michael Paul Wellman.  Spoofing the limit order book: An agent-based model.  In Workshops at the Thirty-First AAAI Conference on Artificial Intelligence, 2017.
  * [14] Jean-Philippe Bouchaud, Marc Mézard, Marc Potters, et al.  Statistical properties of stock order books: empirical results and models.  Quantitative finance, 2(4):251–256, 2002.
  * [15] Matthew Dixon.  Sequence classification of the limit order book using recurrent neural networks.  Journal of computational science, 24:277–286, 2018.
  * [16] Zihao Zhang, Stefan Zohren, and Stephen Roberts.  Deeplob: Deep convolutional neural networks for limit order books.  IEEE Transactions on Signal Processing, 67(11):3001–3012, 2019\.
  * [17] Frédéric Abergel and Aymen Jedidi.  A mathematical approach to order book modeling.  International Journal of Theoretical and Applied Finance, 16(05):1350025, 2013.
  * [18] Jonathan Donier and Julius Bonart.  A million metaorder analysis of market impact on the bitcoin.  Market Microstructure and Liquidity, 1(02):1550008, 2015.
  * [19] Dave Cliff.  Minimal-intelligence agents for bargaining behaviors in market-based environments.  Hewlett-Packard Labs Technical Reports, 1997.


# Thermodynamics of growth in open chemical reaction networks

Shesha Gopal Marehalli Srinivas<EMAIL_ADDRESS>Complex Systems and
Statistical Mechanics, Department of Physics and Materials Science, University
of Luxembourg, L-1511 Luxembourg City, Luxembourg  Francesco Avanzini
<EMAIL_ADDRESS>Complex Systems and Statistical Mechanics,
Department of Physics and Materials Science, University of Luxembourg, L-1511
Luxembourg City, Luxembourg  Department of Chemical Sciences, University of
Padova, Via F. Marzolo, 1, I-35131 Padova, Italy  Massimiliano Esposito
<EMAIL_ADDRESS>Complex Systems and Statistical Mechanics,
Department of Physics and Materials Science, University of Luxembourg, L-1511
Luxembourg City, Luxembourg

###### Abstract

We identify the thermodynamic conditions necessary to observe indefinite
growth in homogeneous open chemical reaction networks (CRNs) satisfying mass
action kinetics. We also characterize the thermodynamic efficiency of growth
by considering the fraction of the chemical work supplied from the
surroundings that is converted into CRN free energy. We find that indefinite
growth cannot arise in CRNs chemostatted by fixing the concentration of some
species at constant values, or in continuous-flow stirred tank reactors.
Indefinite growth requires a constant net influx from the surroundings of at
least one species. In this case, unimolecular CRNs always generate equilibrium
linear growth, i.e., a continuous linear accumulation of species with
equilibrium concentrations and efficiency one. Multimolecular CRNs are
necessary to generate nonequilibrium growth, i.e., the continuous accumulation
of species with nonequilibrium concentrations. Pseudo-unimolecular CRNs – a
subclass of multimolecular CRNs – always generate asymptotic linear growth
with zero efficiency. Our findings demonstrate the importance of the CRN
topology and the chemostatting procedure in determining the dynamics and
thermodynamics of growth.

Significant progress has been achieved over the last decades in establishing a
thermodynamic description of open CRNs both at the stochastic and
deterministic level [1, 2, 3, 4, 5, 6, 7]. The thermodynamic framework
complements the dynamical one by enabling energetic considerations of various
complex chemical phenomena such as the cost of sustaining nonequilibrium
steady states [8], chemical oscillations [9], chaotic dynamics [10], patterns
[11, 12, 13], waves [14, 15]. It has also been used to study the efficiency of
energy transduction [16, 17, 18], energy storage [19] and various forms of
information processing [20, 21]. But surprisingly little attention has been
given to the phenomenon of chemical growth. Under which conditions does an
open CRN show growth by extracting matter and energy from its surroundings and
with what efficiency? The relevance of this question for biology is obvious.
Previous works on growth have mainly focused on describing the growth dynamics
of CRNs with irreversible reactions [22, 23, 24, 25, 26], implicitly assuming
infinite resources to grow. Thus, exponential or even hyperbolic growth was
observed. In the case of reversible CRNs, numerical studies of transient
growth regimes have been considered [27]. In this Letter, we focus on the
question of asymptotic concentration growth by considering the simplest setup
of homogeneous, open, and reversible CRNs in ideal dilute solutions. The
dynamics of such CRNs is given by mass action kinetics. We consider three
different chemostatting mechanisms for opening CRNs: holding the concentration
of some species constant (concentration control), imposing a constant influx
or outflux of some species (flux control), and a constant influx and an
outflux of species proportional to their concentration (mixed control). Many
standard experimental setups can be described in terms of those chemostatting
mechanisms, such as batch control [28], continuous-flow stirred tank
reactors[29], and serial transfer experiments [30]. By analytically studying
growing unimolecular and pseudo-unimolecular CRNs and numerically studying
four multimolecular CRNs, we clarify the important role of the topology of the
CRN and the chemostatting procedure in determining when growth with
equilibrium or nonequilibrium concentrations is possible.

We start by introducing setup and dynamics of open CRNs. Chemical species
$\alpha$ are interconverted by (elementary[31], mass balanced) chemical
reactions $\rho$,

$\bm{\alpha}\cdot\bm{\nu}_{+\rho}\xleftrightharpoons[-\rho]{+\rho}\bm{\alpha}\cdot\bm{\nu}_{-\rho}\
,\ $ (1)

where the vectors $\bm{\nu}_{\pm\rho}$ specify the number of each species
participating in the reaction $\pm\rho$ and
$\bm{\alpha}=(\dots,\alpha,\dots)^{\intercal}$. The topology of CRNs is
encoded in the stoichiometric matrix $\mathbb{S}$, whose columns
$\mathbb{S}_{\rho}\coloneqq\bm{\nu}_{-\rho}-\bm{\nu}_{+\rho}$, quantify the
net change in the number of the species after the reaction $\rho$. The
concentration vector, $\bm{z}(t)=(\dots,[\alpha](t),\dots)^{\intercal}$, obeys
the rate equation

$\frac{d\bm{z}}{dt}=\mathbb{S}\bm{j}(\bm{z})+\bm{I}(\bm{z})\ .\ $ (2)

The first term in Eq. (2) is the product of the stoichiometric matrix
$\mathbb{S}$ and the reaction current vector
$\bm{j}(\bm{z})=(\dots,j_{\rho}(\bm{z}),\dots)^{\intercal}$. It tracks the
total change in concentrations due to chemical reactions. Each reaction
current is expressed as the difference between a forward and a backward flux,
$j_{\rho}(\bm{z})=j_{+\rho}(\bm{z})-j_{-\rho}(\bm{z})$. From mass action
kinetics,

$j_{\pm\rho}(\bm{z})=k_{\pm\rho}\bm{z}^{\nu_{\pm\rho}}\ ,\ $ (3)

where $k_{\pm\rho}$ are rate constants and the notation
$\bm{a}^{\bm{b}}\coloneqq\Pi_{j}a_{j}^{b_{j}}$ is used. The second term
$\bm{I}(\bm{z})$ in Eq. (2) results from chemostatting a subset of species
(denoted Y) and quantifies the exchange currents with the surroundings. The
remaining species (denoted X) are called internal species. As a result, Eq.
(2) splits into $d_{t}\bm{x}=\mathbb{S}^{X}\bm{j}(\bm{z})$ and
$d_{t}\bm{y}=\mathbb{S}^{Y}\bm{j}(\bm{z})+\bm{I}^{Y}(\bm{z})$, where $\bm{x}$
and $\bm{y}$ represent the concentrations of the internal and chemostatted
species, respectively. The three types of chemostatting procedures we consider
are formalized as follows. Concentration control maintains $\bm{y}$ constant
in time: $\bm{I}^{Y}(\bm{z})=-\mathbb{S}^{Y}\bm{j}(\bm{z})$. Flux control
imposes a constant flux of species into or out of the system:
$\bm{I}^{Y}(z)=\tilde{\bm{I}}$. Mixed control imposes a constant influx of
species into the system and extracts species from the system proportionally to
their concentrations:
$\bm{I}^{Y}(z)=-\tilde{\mathbb{D}}\bm{y}+\tilde{\bm{I}}$, where
$\tilde{\mathbb{D}}$ is a diagonal matrix and $\tilde{\bm{I}}$ is a constant
(positive) vector. Continuous-flow stirred tank reactors correspond to a form
of mixed control in which all the species are extracted with the same rate,
i.e., $\tilde{\mathbb{D}}={k}^{e}\mathbb{1}$ [30]. From Eqs. (2) and (3), we
note that unimolecular CRNs, which consist exclusively of reactions of the
form $\alpha\xrightleftharpoons{}\beta$, will follow a linear dynamics for any
chemostatting procedure. CRNs with multimolecular reactions will always give
rise to nonlinear dynamics for flux and mixed control. However, under
concentration control, some of them may give rise to linear dynamics. We call
these pseudo-unimolecular CRNs.

We now turn to the thermodynamics of open CRNs [4, 6, 7]. The chemical
potential for species in homogeneous dilute solutions at temperature $T$ is
given by the vector $\bm{\mu}(\bm{z})=\bm{\mu}^{0}+RT\text{ln}(\bm{z})$, with
$R$ the gas constant and $\bm{\mu}^{0}$ the vector of standard chemical
potentials. From the form of the chemical potential, we deduce the
nonequilibrium Gibbs free energy,
$G(\bm{z})=\sum_{\alpha}\mu_{\alpha}[\alpha]-RT\sum_{\alpha}[\alpha]$, and the
entropy production rate (EPR), $T\dot{\Sigma}=-\bm{\mu}\mathbb{S}\bm{j}$. The
thermodynamics is linked to the dynamics via the local detailed balance,

$RT\leavevmode\nobreak\
\text{ln}\left({k_{+\rho}}/{k_{-\rho}}\right)=-\bm{\mu}^{0}\cdot\mathbb{S}_{\rho}\
.\ $ (4)

This allows us to express the EPR as
$\dot{\Sigma}={R}\sum_{\rho}(j_{+\rho}-j_{-\rho})\text{ln}(j_{+\rho}/j_{-\rho})\geq
0$, and derive a nonequilibrium second law
$T\dot{\Sigma}=\dot{w}_{\text{c}}-d_{t}{G}\geq 0$. The term
$\dot{w}_{\text{c}}=\sum_{\alpha}\mu_{\alpha}I_{\alpha}$, called the chemical
work rate, is the total work done on the CRN per unit time. Thus, the second
law states that the dissipation is the difference between the total work done
on the CRN and the change in its free energy. We further decompose the
chemical work rate $\dot{w}_{c}$ using conservation laws, which are linearly
independent left null eigenvectors $\ell^{\lambda}$ of the stoichiometric
matrix $\mathbb{S}$. They identify moieties, that are parts of (or entire)
molecules that are left intact by the reactions. Their concentrations
${L}^{\lambda}=\ell^{\lambda}\cdot\bm{z}$ are conserved in closed CRNs. In
open CRNs,

$d_{t}{L}^{\lambda}=\ell^{\lambda}\cdot
d_{t}\bm{z}=\ell^{\lambda}\cdot\bm{I}(z)\\\ .$ (5)

If $\ell^{\lambda}\cdot\bm{I}(z)\neq 0$ (resp. $=0$), the moiety concentration
${L}^{\lambda}$ is no longer conserved (resp. still conserved) and the
conservation law $\ell^{\lambda}$ is said to be broken (resp. unbroken). Note
that every CRN has at least one conservation law that involves all species,
denoted $\ell^{m}$, with the corresponding moiety concentration
$L^{m}=\ell^{m}\cdot\bm{z}$ called the mass density [32], which is always
broken in an open CRN.

Since the moiety concentrations ${L}^{\lambda}$ change only due to the
exchange currents, we split the chemical work rate into the work done in
changing ${L}^{\lambda}$ and the remaining. However, as chemostatting a
species does not always break a conservation law, we first split the set of
chemostatted species $Y$ into potential species $Y_{p}$ which break
conservation laws and force species $Y_{f}=Y\setminus Y_{p}$. This allows us
to associate a moiety to a single $Y_{p}$ species and express its
concentration according to
$\bm{m}(\bm{z})=\big{(}\mathbb{L}^{b}_{Y_{p}}\big{)}^{-1}\mathbb{L}^{b}\cdot\bm{z}$,
where $\mathbb{L}^{b}$ is a matrix with broken conservation laws as rows while
$\mathbb{L}^{b}_{Y_{p}}$ is its square submatrix [6]. The work done in
changing the exchanged moiety concentrations, called the moiety work rate,
then reads $\dot{w}_{\text{m}}=\bm{\mu}_{Y_{p}}d_{t}\bm{m}$. The remaining
part of the chemical work rate is called the nonconservative work rate,
$\dot{w}_{\text{nc}}=\dot{w}_{\text{c}}-\dot{w}_{\text{m}}$. The
nonconservative work rate identifies the energetic cost of maintaining fluxes
of the same moiety between the various chemostats [4, 6] and can be rewritten
as
$\dot{w}_{\text{nc}}=\big{(}\bm{\mu}_{Y}-\bm{\mu}_{Y_{p}}\big{(}\mathbb{L}^{b}_{Y_{p}}\big{)}^{-1}\mathbb{L}^{b}\big{)}\cdot\bm{I}^{Y}$.
As a result, the nonequilibrium second law of CRNs can be rewritten in a form
particularly well suited to analyzing growth:

$T\dot{\Sigma}=\dot{w}_{\text{nc}}+\dot{w}_{\text{m}}-d_{t}{G}\geq 0\;.$ (6)

Closed CRNs always reach a fixed point, $\bm{z}_{\text{eq}}$, called the
equilibrium state, such that $\bm{j}(\bm{z}_{\text{eq}})=0$ and
$T\dot{\Sigma}=\dot{w}_{\text{nc}}=0$. Fixed points in open CRNs can either be
equilibria or nonequilibrium steady states (NESS), $\bm{z}_{\text{ss}}$, with
$\bm{j}(\bm{z}_{\text{ss}})\neq 0$ and $T\dot{\Sigma}=\dot{w}_{\text{nc}}>0$.
Note, for any fixed point of a CRN, $\dot{w}_{\text{m}}=d_{t}G=0$. Often,
chemostatting CRNs is done by concentration control. Generically, such open
CRNs relax to a NESS or exhibit complex behaviors like limit cycles or
multistability. However, no asymptotic growth was observed till recently, in
Ref. [7], where numerically, under both flux and mixed control, a CRN was
shown to tend towards a state with infinite concentrations. In this work, we
investigate the dynamics and energetics of growing CRNs.

Indefinite growth implies that at least one species keeps accumulating in the
system. It is thus equivalent to the statement that
$\lim_{t\to\infty}||\bm{z}(t)-\bm{z}(0)||=\infty$ (where $||\bullet||$ denotes
the Euclidean norm), or that the mass density diverges
$\lim_{t\to\infty}L^{m}(t)=\infty$. Thermodynamically, the accumulation of
species results in free energy that diverges in time. This tells us that the
work needed from the surroundings also has to diverge which immediately
implies that closed CRNs do not grow and only open CRNs may grow. The
chemostatting procedure plays a key role in determining the growth conditions.
For example, in the case of continuous-flow stirred tank reactors, Eq. (5) for
the mass density becomes
$d_{t}{L}^{m}=-{k}^{e}{L}^{m}+\ell^{m}_{Y}\cdot\tilde{\bm{I}}$, implying a
bounded solution and therefore no growth for any choice of rate constants
$\\{k_{\pm\rho},{k}^{e}\\}$ and current $\tilde{\bm{I}}$ [30]. As another
example, under flux control, all CRNs show growth if
$d_{t}{L}^{m}=\ell^{m}\cdot\tilde{\bm{I}}>0$, with all growing concentrations
diverging at most linearly in time.

For a systematic analysis, we first consider CRNs that are dynamically linear
and can thus be characterized analytically. We note that the detailed
calculations and a more extensive analysis are reported in a companion paper
[33]. The first important result is that unimolecular CRNs can produce
equilibrium growth under flux control and no growth otherwise. For any
chemostatting procedure of a unimolecular CRN, Eq. (2) can be rewritten in the
form

$d_{t}\bm{a}=\mathbb{W}\bm{a}-\mathbb{D}\bm{a}+\bar{\bm{I}}\,.$ (7)

Here, $\bm{a}=\bm{z}$ under flux and mixed control, while $\bm{a}=\bm{x}$
under concentration control. $\mathbb{W}$ is a rate matrix [34] whose
nondiagonal entries are given by
$\mathbb{W}_{\alpha,\beta}=\\{\sum_{\rho}k_{+\rho}\leavevmode\nobreak\
\leavevmode\nobreak\ \big{|}\leavevmode\nobreak\
\rho:\beta\xrightleftharpoons{+\rho}\alpha\\}$ and diagonal entries are such
that the sum over each column is zero. Due to the local detailed balance
condition (4), $\mathbb{W}$ is a detailed balanced matrix with a zero
eigenvector proportional to $\exp(-{\bm{\mu^{0}}}/{{RT}})$. Under all kinds of
control, $\bar{\bm{I}}$ is a constant vector. Under mixed control or
concentration control, $\mathbb{D}$ is a non-negative diagonal matrix with at
least one nonzero element. As a result, all the eigenvalues of the matrix
$\mathbb{W}-\mathbb{D}$ are negative and growth cannot occur.

Flux control instead corresponds to $\mathbb{D}=0$ in Eq. (7). The mass
density grows linearly in time at a rate given by the total injection current,
$d_{t}L^{m}(t)=\sum_{\alpha}\bar{{I}}_{\alpha}$, and the growth of all species
is controlled by it at large times:

$[\alpha](t)=\underbrace{\left(\frac{e^{{-\mu^{0}_{\alpha}}/{{RT}}}}{\sum_{\alpha^{\prime}}e^{{-\mu^{0}_{\alpha^{\prime}}}/{{RT}}}}\right)L^{m}(t)}_{=[\alpha]_{\text{eq}}\left(L^{m}(t)\right)}+\leavevmode\nobreak\
c_{\alpha}(\bar{\bm{I}})=\mathcal{O}(t)\ .\ $ (8)

Here, $\\{c_{\alpha}(\bar{\bm{I}})\\}$ are time-independent coefficients and
linear functions of the injection currents, while
$[\alpha]_{\text{eq}}\left(L^{m}(t)\right)$ is the equilibrium value to which
the CRN would relax if it were closed (by stopping the injection) at time $t$.
This means that the growth process follows, up to the offset
$c_{\alpha}(\bar{\bm{I}})$, that equilibrium value. Turning to the second law,
at large times one finds that

$\displaystyle T\dot{\Sigma}$
$\displaystyle\sim\dot{w}_{\text{nc}}=\mathcal{O}\left(t^{-1}\right),$ (9)
$\displaystyle d_{t}{G}$
$\displaystyle\sim\dot{w}_{\text{m}}=\mathcal{O}\left(\text{ln}\left(t\right)\right).$
(10)

The physical interpretation of this result is simple. Since the CRN
concentrations scale extensively over time, while the offset in Eq. (8)
remains intensive, growth over large times resembles a quasistatic process
where the total dissipation becomes negligible and the moeity work performed
by the injection is reversibly converted into free energy. Defining the
thermodynamic efficiency with which the chemical work is converted into free
energy, $\eta=d_{t}G/\dot{w}_{\text{c}}$, and using Eqs. (10), (9) and (6), we
see that $\eta$ tends to one at large times for unimolecular CRNs, indicating
reversible energy conversions. In other words, equilibrium growth can be seen
at large times as a moving equilibrium following the mass density displaced by
the injection.

We now analyze pseudo-unimolecular CRNs. These are multimolecular CRNs where
the concentrations of certain species, called hidden species, are held
constant via concentration control, and where the dynamics of the remaining
concentrations, called dynamically linear species, become linear. For example,
consider the single reaction
$\text{E}+\text{S}\xrightleftharpoons{}\text{ES}$. If the concentration
$[\text{S}]$ is held constant, the dynamics becomes linear. By absorbing the
(constant) concentrations of the hidden species into the rate constants, the
rate equation for pseudo-unimolecular CRNs can still be written as Eq. (7)
with a rate matrix $\mathbb{W}$ that is not detailed balanced. For the same
reasons as for unimolecular CRNs, mixed control and concentration control on
the dynamically linear species will not give rise to any growth. Instead, flux
control on the dynamically linear species will produce growth which over long
times reads

$[\alpha](t)=\underbrace{\pi^{\text{ss}}_{\alpha}L^{m}(t)}_{=[\alpha]_{\text{ss}}\left(L^{m}(t)\right)}+\leavevmode\nobreak\
c_{\alpha}.$ (11)

Here, $\pi^{\text{ss}}_{\alpha}$ is a component of the normalized right null
vector of $\mathbb{W}$ and $[\alpha]_{\text{ss}}\left(L^{m}(t)\right)$ is the
steady state (generically a NESS) to which the CRN would relax if it were
closed (by stopping the injection) at time $t$. In the long time limit, the
thermodynamic quantities are of the form,

$\displaystyle T\dot{\Sigma}$
$\displaystyle\sim\dot{w}_{\text{nc}}=\mathcal{O}(t)\,,$ (12) $\displaystyle
d_{t}{G}$
$\displaystyle\sim\dot{w}_{\text{m}}=\mathcal{O}\left(\text{ln}\left(t\right)\right)\;.$
(13)

Exactly as for growing unimolecular CRNs (see Eq. (10)) the moiety work rate
performed by the injection is converted into free energy. The key difference
is that this conversion is generically irreversible, and thus gives rise to an
extensive EPR in time. In contrast to the unimolecular case, the efficiency
$\eta=d_{t}G/\dot{w}_{\text{c}}$ tends to zero at large times (using Eqs.
(12), (13) and (6)) indicative of very inefficient irreversible conversions.
This non-equilibrium growth process can be seen over long times as a moving
NESS following the mass density displaced by the injection. The dissipation is
essentially the nonconservative work dissipated in the NESS due to the
presence of the chemostatted hidden species.

Figure 1: Evolution of a) the concentrations, b) the contributions to the EPR
in Eq. (6) and the efficiency $\eta$, for the Michaelis Menten CRN in a) under
mixed control when S is injected at the rate $I_{\ch S}=1$ while P is
extracted at the rate $-{k}^{e}_{\ch{P}}[\ch{P}]$. We rescaled time,
concentrations, and energy/time by $1/k_{-1}$, $k_{-1}/k_{+1}$ and
$RTk_{-1}^{2}/k_{+1}$, respectively. $[\text{E}](0)=0.3$,
$[\text{ES}](0)=0.1$, $[\text{S}](0)=1.9$, $[\text{P}](0)=0.4$,
${k}^{e}_{\ch{P}}=0.5$, $k_{\pm 1}=k_{\pm 2}=1$, $\mu^{0}_{\ch{E}}=1$,
$\mu^{0}_{\ch{ES}}=2$, $\mu^{0}_{\ch{S}}=1$ and $\mu^{0}_{\ch{P}}=1$.

In summary, by considering dynamically linear CRNs, we prove that only flux
chemostatting generates growth in such CRNs. Unimolecular CRNs only display
equilibrium growth, while pseudo-unimolecular CRNs can give rise to
nonequilibrium growth. A striking feature of these CRNs is a splitting of the
second law, where the total dissipation is essentially the nonconservative
work rate whereas the moiety work rate is essentially converted into free
energy.

We now turn to multimolecular CRNs. General statements are notoriously
difficult to make in this case. However, it is conjectured that multimolecular
CRNs cannot grow under concentration control. Although this is proved for the
cases of single linkage class CRNs [35] and strongly endotactic CRNs [36], a
general proof is still lacking. We considered four different multimolecular
CRNs, based on their different topological properties, which are analyzed
(analytically and numerically) in detail in the companion paper [33]. Two of
them are shown in the present letter. All support the conjecture. The
intuition behind this conjecture is that the accumulation of any species will
eventually produce reaction currents (due to mass action kinetics) which will
consume the species and these currents would be sufficiently large to balance
the external currents generated by the chemostats. Using these four models, we
also found that, contrary to dynamically linear CRNs, some multimolecular CRNS
can grow under mixed control. Indeed, in Fig. 1a, we show that the Michaelis
Menten CRN under mixed control grows out of equilibrium.

Figure 2: Evolution of a) the concentrations, b) the contributions to the EPR
in Eq. (6) and the efficiency $\eta$, for the autocatalytic CRN in a) under
flux control when S is injected at the rate $I_{\ch S}=1$ while P is extracted
at the rate $I_{\ch P}=-0.1$. Here, we rescale time, concentrations, and
energy/time by $1/k_{+1}$, $k_{+1}/k_{+3}$ and $RTk_{+1}^{2}/k_{+3}$,
respectively. $[\ch{E}](0)=1.3$, $[\ch{E}^{*}](0)=0.1$, $[\ch{F}](0)=1.9$,
$[\ch{W}](0)=0.4$, $k_{\pm 1}=k_{\pm 2}=k_{\pm 3}=1$, $\mu^{0}_{\ch{E}}=1$,
$\mu^{0}_{\ch{E^{*}}}=2$, $\mu^{0}_{\ch{S}}=1$, $\mu^{0}_{\ch{P}}=2$.

For all four models, irrespective of the type of control, in growth regimes,
we observed that there exists a choice of potential species such that the
thermodynamics displays the same splitting of the second law seen in
unimolecular and pseudo-unimolecular CRNs (see Eqs. (9), (10), (12) and (13))
of the form

$\displaystyle T\dot{\Sigma}$ $\displaystyle\sim\dot{w}_{\text{nc}}\geq 0\,,$
(14) $\displaystyle d_{t}G$ $\displaystyle\sim\dot{w}_{\text{m}}>0\,.$ (15)

However, the generality of this observation remains to be tested beyond these
models. Furthermore, the energy conversions in growing multimolecular CRNs are
generically such that the efficiency might be strictly between zero and one
(see inset of Fig. 1b). Finally, multimolecular CRNs can display growth
regimes that are significantly different from those in dynamically linear
CRNs. Only a subset of species may grow and the temporal scalings of the
concentrations may differ, as shown in Fig. 2a.

Figure 3: Long-time behavior of concentrations, EPR, and efficiency. No growth
under orange (resp. yellow) conditions (resp., but definite proof is still
missing). Equilibrium (resp. nonequilibrium) growth under solid (resp.
striped) blue conditions.

Conclusion and perspectives \- Our main findings are summarized in Fig. 3. We
considered growth in dilute ideal solutions, implicitly assuming that over the
timescales considered, the species concentrations remain negligible compared
to the solvent. Approaches to growth introducing an equation of state coupling
volume and species concentrations have been considered but only for detailed
balanced dynamics (which always relax to equilibrium) [37, 38, 39]. A future
perspective may be to consider the energetics of such systems for non-detailed
balanced dynamics. Characterizing the effect of diffusion or the effect of
compartment division [40, 41, 42] on growth may also be considered. Finally,
we observed that autocatalysis does not seem to play an important role on the
long time growth (see companion paper [33]). Perhaps the role of autocatalysis
is more apparent in the short-time growing behavior [43, 44, 27].

Acknowledgements \- The authors would like to thank Danilo Forastiere and
Gianmaria Falasco for helpful discussions. This research was supported by the
Luxembourg National Research Fund (FNR), via the research funding schemes
PRIDE (Grant No. 19/14063202/ACTIVE), CORE project ChemComplex (Grant No.
C21/MS/16356329), and by project INTER/FNRS/20/15074473 funded by F.R.S.-FNRS
(Belgium) and FNR (Luxembourg).

## References

  * Qian and Beard [2005] H. Qian and D. A. Beard, Biophysical Chemistry 114, 213 (2005).
  * Schmiedl and Seifert [2007] T. Schmiedl and U. Seifert, J. Chem. Phys. 126, 044101 (2007).
  * Polettini and Esposito [2014] M. Polettini and M. Esposito, J. Chem. Phys. 141, 024117 (2014).
  * Rao and Esposito [2016] R. Rao and M. Esposito, Phys. Rev. X 6, 041064 (2016).
  * Rao and Esposito [2018] R. Rao and M. Esposito, J. Chem. Phys. 149, 245101 (2018).
  * Avanzini  _et al._ [2021] F. Avanzini, E. Penocchio, G. Falasco, and M. Esposito, J. Chem. Phys. 154, 094114 (2021).
  * Avanzini and Esposito [2022] F. Avanzini and M. Esposito, J. Chem. Phys. 156, 014116 (2022).
  * Gaspard [2004] P. Gaspard, J. Chem. Phys. 120, 8898 (2004).
  * Andrieux and Gaspard [2008] D. Andrieux and P. Gaspard, J. Chem. Phys. 128, 154506 (2008).
  * Gaspard [2020] P. Gaspard, Chaos 30, 113103 (2020).
  * Falasco  _et al._ [2018] G. Falasco, R. Rao, and M. Esposito, Phys. Rev. Lett. 121, 108301 (2018).
  * Avanzini  _et al._ [2020] F. Avanzini, G. Falasco, and M. Esposito, Phys. Rev. E 101, 060102 (2020).
  * Aslyamov  _et al._ [2023] T. Aslyamov, F. Avanzini, E. Fodor, and M. Esposito, Phys. Rev. Lett. 131, 138301 (2023).
  * Avanzini  _et al._ [2019] F. Avanzini, G. Falasco, and M. Esposito, J. Chem. Phys. 151, 234103 (2019).
  * Kumar and Gangopadhyay [2021] P. Kumar and G. Gangopadhyay, Phys. Rev. E 104, 014221 (2021).
  * Hill [1966] T. L. Hill, J. Theor. Biol. 10, 442 (1966).
  * Wachtel  _et al._ [2022] A. Wachtel, R. Rao, and M. Esposito, J. Chem. Phys. 157, 024109 (2022).
  * Corra  _et al._ [2022] S. Corra, M. T. Bakić, J. Groppi, M. Baroncini, S. Silvi, E. Penocchio, M. Esposito, and A. Credi, Nat. Nanotech. 17, 746 (2022).
  * Penocchio  _et al._ [2019] E. Penocchio, R. Rao, and M. Esposito, Nat. Commun. 10, 3865 (2019).
  * Amano  _et al._ [2022] S. Amano, M. Esposito, E. Kreidt, D. A. Leigh, E. Penocchio, and B. M. W. Roberts, Nat. Chem. 14, 530 (2022).
  * Penocchio  _et al._ [2022] E. Penocchio, F. Avanzini, and M. Esposito, J. Chem. Phys. 157, 034110 (2022).
  * Lin  _et al._ [2020] W.-H. Lin, E. Kussell, L.-S. Young, and C. Jacobs-Wagner, PNAS 117, 27795 (2020).
  * Szathmáry [1991] E. Szathmáry, Trends Ecol. Evol. 6, 366 (1991).
  * Wills  _et al._ [1998] P. R. Wills, S. A. Kauffman, B. M. Stadler, and P. F. Stadler, Bull. Math. Biol. 60, 1073 (1998).
  * Iyer-Biswas  _et al._ [2014] S. Iyer-Biswas, G. E. Crooks, N. F. Scherer, and A. R. Dinner, Phys. Rev. Lett. 113, 028101 (2014).
  * Nandori and Young [2022] P. Nandori and L.-S. Young, PNAS 119, e2214282119 (2022).
  * Sarkar and England [2019] S. Sarkar and J. L. England, Phys. Rev. E 100, 022414 (2019).
  * Sorrenti  _et al._ [2017] A. Sorrenti, J. Leira-Iglesias, A. Sato, and T. M. Hermans, Nat. Commun. 8, 15899 (2017).
  * Liu  _et al._ [2023] K. Liu, A. Blokhuis, C. van Ewijk, A. Kiani, J. Wu, W. H. Roos, and S. Otto, Nat. Chem. (2023), 10.1038/s41557-023-01301-2.
  * Blokhuis  _et al._ [2018] A. Blokhuis, D. Lacoste, and P. Gaspard, J. Chem. Phys. 148, 144902 (2018).
  * Svehla [1993] G. Svehla, Pure and applied chemistry 65, 2291 (1993).
  * Angeli [2009] D. Angeli, in  _2009 European Control Conference (ECC)_ (2009) pp. 649–657.
  * Srinivas  _et al._ [2023] S. M. Srinivas, F. Avanzini, and M. Esposito, “Characterizing dynamics and thermodynamics of growth in open chemical reaction networks,” (2023), submitted to PRE.
  * Van Kampen [2007] N. G. Van Kampen,  _Stochastic processes in physics and chemistry_ (North Holland, 2007).
  * Anderson [2011] D. F. Anderson, J. Math. Chem. 49, 2275 (2011).
  * Gopalkrishnan  _et al._ [2014] M. Gopalkrishnan, E. Miller, and A. Shiu, SIAM J. Appl. Dyn. 13, 758 (2014).
  * Sughiyama  _et al._ [2022] Y. Sughiyama, A. Kamimura, D. Loutchko, and T. J. Kobayashi, Phys. Rev. Res. 4, 033191 (2022).
  * Bigan  _et al._ [2015] E. Bigan, J. Steyaert, and S. Douady, Electronic Notes in Theoretical Computer Science 316, 3 (2015), 5th International Workshop on Static Analysis and Systems Biology (SASB 2014).
  * Kondo and Kaneko [2011] Y. Kondo and K. Kaneko, Phys. Rev. E 84, 011927 (2011).
  * Bauermann  _et al._ [2022] J. Bauermann, C. A. Weber, and F. Jülicher, Ann. Phys. 534, 2200132 (2022).
  * Klosin  _et al._ [2020] A. Klosin, F. Oltsch, T. Harmon, A. Honigmann, F. Jülicher, A. A. Hyman, and C. Zechner, Science 367, 464 (2020).
  * Ianeselli  _et al._ [2022] A. Ianeselli, D. Tetiker, J. Stein, A. Kühnlein, C. B. Mast, D. Braun, and T.-Y. Dora Tang, Nat. Chem. 14, 32 (2022).
  * Unterberger and Nghe [2022] J. Unterberger and P. Nghe, J. Math. Biol. 85, 26 (2022).
  * Blokhuis  _et al._ [2020] A. Blokhuis, D. Lacoste, and P. Nghe, PNAS 117, 25230 (2020).