# A Simple Comparison of Biochemical Systems Theory and Metabolic Control
Analysis

Herbert M. Sauro
Department of Bioengineering
University of Washington
Seattle, WA 98195-5061
<EMAIL_ADDRESS>
###### Abstract

This paper explores some basic concepts of Biochemical Systems Theory (BST)
and Metabolic Control Analysis (MCA), two frameworks developed to understand
the behavior of biochemical networks. Initially introduced by Savageau, BST
focuses on system stability and employs power laws in modeling biochemical
systems. On the other hand, MCA, pioneered by authors such as Kacser and Burns
and Heinrich and Rapoport, emphasizes linearization of the governing equations
and describes relationships (known as theorems) between different measures.
Despite apparent differences, both frameworks are shown to be equivalent in
many respects. Through a simple example of a linear chain, the paper
demonstrates how BST and MCA yield identical results when analyzing steady-
state behavior and logarithmic gains within biochemical pathways. This
comparative analysis highlights the interchangeability of concepts such as
kinetic orders, elasticities and other logarithmic gains.

_Keywords_ Biochemical Systems Theory $\cdot$ Metabolic Control Analysis
$\cdot$ Logarithmic gain

##  1 Introduction

In the late 60s and early 1970s Michael Savageau (Savageau, 1972, 1976)
developed an approach called biochemical systems theory (BST). This was a set
of mathematical tools to help understand the behavior of biochemical networks.
This work arose alongside another similar and independently developed approach
in Europe called metabolic control analysis (MCA) by authors Kacser and Burns
(Kacser, 1973) and Heinrich and Rapoport (Heinrich and Rapoport, 1974). Both
approaches appear at first glance to be different but are actually identical
in many respects. Some of the differences arise from different notation but
there are also differences in emphasis. BST focuses more on system stability
and the use of power laws, while MCA barely mentions stability and is based
instead on a direct linearization of the governing equations. Both approaches
examine how behavior at the systems level is governed by properties at the
local level. MCA also emphasizes theorems that describe relationships between
the different measures. Despite their differences, both approaches have
provided a wealth of insight into how biochemical networks operate with BST
providing a fairly thorough analysis of systems with negative feedback. In
this short article, I will show, using some simple examples, the equivalence
between the two approaches. No significant new results are presented and the
article serves more as an introduction to these two important frameworks. We
will start by considering BST.

A key innovation with BST was the use of power laws as reasonable
approximations to use when building models of biochemical networks. For
example, given the following single irreversible reaction with substrate
$X_{1}$:

$X_{1}\stackrel{{\scriptstyle v_{i}}}{{\longrightarrow}}$

the reaction rate, $v_{i}$, can be written in a power law form:

$v_{i}=\alpha_{i}X_{1}^{g_{i1}}$

where $\alpha_{i}$ is called the apparent rate constant and the exponent
$g_{i1}$ is called the apparent kinetic order. If the kinetic order is one,
the equation reduces to the simple first-order mass-action rate law. The power
law formulation arises from a Taylor expansion in log space (Savageau, 1969)
at a specific operating in the kinetic response. The kinetic orders are local
sensitivities in log space, that is:

$g_{ij}=\frac{\partial v_{i}}{\partial
X_{j}}\frac{X_{j}}{v_{i}}=\frac{\partial\log(v_{i})}{\partial\log(X_{j})}$

These are exactly equivalent to the elasticities used in MCA.

BST also uses the symbol $h$ for kinetic orders specifically for degradation
steps however this distinction is less important as we’ll see when building
complete models. Given this formalism is it straightforward to write a set of
differential equations for a given model where we assign each reaction a
specific power law. For example, consider a four-step linear chain:

$X_{o}\stackrel{{\scriptstyle
v_{1}}}{{\longrightarrow}}X_{1}\stackrel{{\scriptstyle
v_{2}}}{{\longrightarrow}}X_{2}\stackrel{{\scriptstyle
v3}}{{\longrightarrow}}X_{3}\stackrel{{\scriptstyle
v_{4}}}{{\longrightarrow}}$

We will assume that the species $X_{o}$ is fixed so that the pathway can
sustain a steady-state. The differential equations for this system can be
easily written in terms of the production and consumption rates as:

$\displaystyle\begin{split}\frac{dX_{1}}{dt}&=v_{1}-v_{2}\\\\[4.0pt]
\frac{dX_{2}}{dt}&=v_{2}-v_{3}\\\\[4.0pt]
\frac{dX_{3}}{dt}&=v_{3}-v_{4}\end{split}$ (1)

For each $v_{i}$ we can substitute the appropriate power law to give:

$\displaystyle\begin{split}\frac{dX_{1}}{dt}&=\alpha_{1}X_{o}^{g_{10}}-\beta_{1}X_{1}^{h_{11}}\\\\[4.0pt]
\frac{dX_{2}}{dt}&=\beta_{1}X_{1}^{h_{11}}-\beta_{2}X_{2}^{h_{22}}\\\\[4.0pt]
\frac{dX_{3}}{dt}&=\beta_{2}X_{2}^{h_{22}}-\beta_{3}X_{2}^{h_{33}}\end{split}$
(2)

The sign of the kinetic order depends on whether the effector increases
(positive) or decreases (negative) a reaction rate. This means that the
kinetic orders for inhibitors will be negative in value. In this example they
are all positive because increases in substrate concentration will increase
the reaction rate.

Note the use of $h$ as the kinetic order for the degradation rate for $v_{2}$.
That is, the symbol $h$ isn’t always used to describe degradation steps but,
as in this case, it is also used in a consumption step when describing the
rate of change of $X_{2}$. The second equation, $dX_{2}/dt$, has the same $h$
kinetic order because the same reaction is now a production rate for $X_{2}$.
In fact, in all these equations, the $g$ kinetic orders, other than $g_{10}$
are replaced by $h$ type kinetic orders irrespective of whether we are dealing
with a production or consumption step. In general, the distinction between
production and consumption kinetic orders does not appear to be so important,
especially for reaction networks where each reaction step is explicitly
written in the differential equation. Additionally, the notation changes
depending on what formalism is being used since BST has a variety of
representations it can employ. The two main approaches that are used are the
Generalized Mass-Action (GMA) and the S-Systems (Savageau, 1988) formalism.
When dealing with linear chains, as we will do here, the two approaches are
identical. They diverge when a model has branches or cycles. In these
situations, GMA has explicit terms for each reaction rate. For example, for a
simple branch, Figure 1, the rate of change of the branch species will
comprise three power laws, one for the production rate and two for the
consumption rates emerging from the branch. Equation (3) shows a GMA
representation of the rate of change of the branch species $X_{1}$. For GMA
models, the kinetics orders are generally represented using the $f$ symbol. A
more complex example that uses the GMA approach, can be found at Alves et al.
(2008).

Figure 1: Simple branch pathway with two rates $v_{2}$ and $v_{3}$ consuming
$X_{1}$.

$\displaystyle\frac{dX_{1}}{dt}=\alpha_{1}X_{0}^{g_{10}}-\beta_{1}X_{1}^{h_{21}}-\beta_{2}X_{1}^{h_{31}}$
(3)

For an S-System, which is where most of the utility of BST originates, all
power laws associated with the production of a particular species are
aggregated. Likewise, all consumption rates are aggregated. This is where $g$
is more likely to be used for the aggregated production rate and $h$ for the
aggregated consumption steps to keep track of what is what. Aggregation,
however, can be challenging for complex systems that involve branches and
cycles. This may involve, for example, computing the weighted kinetic orders
(Voit, 2013). For linear chains, the S-system formalism is easy to apply since
no aggregation is necessary. This is why GMA and S-Systems are equivalent for
linear chains.

##  2 Examples

To illustrate the utility of BST and how it matches MCA, three examples will
be considered. The first example is the four-step linear chain of reactions we
considered previously. To keep things simple, we will assume that each
reaction is irreversible. That is, the product has no influence over its own
production. This will mean that each power law only has a single species in
the equation which corresponds to the substrate of the reaction. This makes
the analysis a little easier. To make a reaction sensitive to product we only
need to augment the power law with a product term. For example, to make the
first reaction sensitive to product $X_{1}$, we can add a $X_{1}$ term with
its own kinetic order, $g_{11}$:

$v_{1}=\alpha_{1}X_{o}^{g_{10}}X_{1}^{g_{11}}$

We will look at this situation in the next section. Making reactions fully
reversible, in the sense that the reaction rate can go negative, is more
complex and requires choices to be made when aggregating (Sorribas and
Savageau, 1989; Voit, 2013) which we will not consider here.

If we assume the pathway is at steady-state we can set the left-hand side of
(2) to zero:

$\displaystyle\begin{split}0&=\alpha_{1}X_{o}^{g_{10}}-\beta_{1}X_{1}^{h_{11}}\\\\[4.0pt]
0&=\beta_{1}X_{1}^{h_{11}}-\beta_{2}X_{2}^{h_{22}}\\\\[4.0pt]
0&=\beta_{2}X_{2}^{h_{22}}-\beta_{3}X_{2}^{h_{33}}\end{split}$ (4)

Rearrangement and taking logarithms on both sides (Savageau, 1969, 1971)
results in a set of linear equations that can be solved for $\log X_{i}$. If
we designate $\log X_{i}$ with the symbol $y_{i}$, these solutions can be
written as follows:

$\displaystyle\begin{split}y_{1}&=(g_{10}y_{0}-b_{1})/h_{11}\\\\[3.0pt]
y_{2}&=(g_{10}y_{0}-b_{1}-b_{2})/h_{22}\\\\[3.0pt]
y_{3}&=(g_{10}y_{0}-b_{1}-b_{2}-b_{3})/h_{22}\end{split}$ (5)

The $b_{i}$ terms are the log of the ratio of $\beta_{i}$ to $\alpha_{i}$ and
$y_{i}$ as indicated before are the log of the $X_{i}$ variables. We note
again that in this example it is assumed that each reaction is not influenced
by its product. The log transformation of the power law representation is a
key innovation in BST and allows one to derive analytically the steady-state
solution from the original differential equations. This formulation is unique
to BST.

The fact that the $y_{i}$ terms are the log of the $X_{i}$, makes it very easy
to compute the logarithmic gains by simple differentiation of these
expressions. Both BST and MCA focus much of their attention on deriving
logarithmic gains and it is worth briefly saying something about such
quantities. Interestingly the term ‘logarithmic gain’ is not widely used in
the literature even though it expresses accurately the concept. The earliest
reference I could find to this term is from one of Savageau’s earliest papers
(Savageau, 1971). Similar concepts to the logarithmic gain appear in economics
in the form of supply and demand elasticities which is also the origin of the
term ‘elasticity’ in MCA. A logarithmic gain is in simple terms a ratio of
relative changes. In differential calculus, we are familiar with the ratio of
absolute changes, which in the limit become differentials. In a similar
fashion logarithmic gains are defined as the limit of infinitesimal relative
changes. Experimentally, logarithmic gains are much more useful because they
are unit-less and experiments in cell and molecular biology are often
expressed as ratios of unit-less fold changes. Given a variable $Y$ and input
variable $X$. the logarithmic gain of how $X$ can influence $Y$ is described
as a logarithmic gain in a variety of forms:

$\text{logarithmic gain }=\frac{dY}{dX}\frac{X}{Y}=\frac{d\log Y}{d\log
X}\approx\frac{Y\%}{X\%}$

In BST, the logarithmic gain is also symbolized using $L_{YX}$ and in MCA as
$L^{Y}_{X}$ where $L$ is a different letter depending on the gain under
consideration. Numbers may also be used instead of letters for both BST and
MCA. BST also uses an additional gain called the sensitivity. This measures
how a given system property such as a concentration is affected by a change in
a system parameter such as a kinetic order. This is symbolized using $S_{Lh}$.
There is not standardized notation but generally, the aforementioned notation
is commonly found in the literature.

In log form, it is easy to derive logarithmic gain if the equation under study
is already in log form. This is the case with the transformation done by BST.

For example, if we wish to investigate how $X_{o}$ influences the steady-state
levels of $X_{i}$, we only have to differentiate each $y_{i}$ with respect to
$y_{o}$, because by design $y_{o}=\log X_{o}$ and $y_{i}=\log X_{i}$. That is:

$\displaystyle\frac{\partial y_{1}}{\partial y_{0}}=\frac{\partial\log
X_{1}}{\partial\log X_{o}}$ (6)

If we differentiate $y_{1}=(g_{10}y_{0}-b_{1})/h_{11}$ with respect to
$y_{0}$, we obtain (noting that the $b$ terms are constant):

$\displaystyle\frac{\partial\log X_{1}}{\partial\log
X_{o}}=\frac{g_{10}}{h_{11}}$ (7)

As mentioned before logarithmic gains are useful because they can be
interpreted as ratios of percentage changes which makes them experimentally
measurable. The other two logarithmic gains can be derived in the same way to
give the complete set of relations between $X_{o}$ and the three
intermediates:

$\displaystyle\begin{split}\frac{\partial\log X_{1}}{\partial\log
X_{o}}&=\frac{g_{10}}{h_{11}}\\\\[4.0pt] \frac{\partial\log
X_{2}}{\partial\log X_{o}}&=\frac{g_{10}}{h_{22}}\\\\[4.0pt]
\frac{\partial\log X_{3}}{\partial\log
X_{o}}&=\frac{g_{10}}{h_{33}}\end{split}$ (8)

From these equations, we see a clear pattern that makes them easy to extend to
any length pathway. These results are well known (Savageau, 1972) and
illustrate one of the key features of BST which is the ability to find
analytical solutions for the steady-state using simple linear algebra from
which the logarithmic gains are easily derived. This approach works extremely
well for linear chains and linear chains with any manner of feedback loops.

Because the kinetic orders are positive values, each of the logarithmic gains
is also positive. This simply means that increasing $X_{0}$ will increase the
steady-state concentrations of $X_{1},X_{2}$, and $X_{3}$ which makes logical
sense.

### Comparison with MCA

We now compare the BST derivation with that from MCA. The first difference is
that in MCA (Kacser et al., 1995; Heinrich and Rapoport, 1974; Sauro, 2018)
the kinetic orders are called the elasticity coefficients and like BST,
logarithmic gains can be derived in terms of the elasticities.

Over the years a variety of ways have been developed to derive the logarithmic
gains (Sauro, 2018). One of these approaches is given in the appendix which
uses implicit differentiation of the differential equations when written in
functional form. The solutions using this method are given below:

$\displaystyle
R^{X_{1}}_{X_{0}}=\frac{\varepsilon^{1}_{0}}{\varepsilon^{2}_{1}},\qquad
R^{X_{2}}_{X_{0}}=\frac{\varepsilon^{1}_{0}}{\varepsilon^{3}_{2}},\qquad
R^{X_{3}}_{X_{0}}=\frac{\varepsilon^{1}_{0}}{\varepsilon^{4}_{3}}$ (9)

We can compare these equations with the ones derived by BST (8). It is clear
there is a direct equivalence. The kinetic order $g_{10}$ is equivalent to the
elasticity $\varepsilon^{1}_{0}$ and the kinetic orders $h_{11},h_{22}$, and
$h_{33}$ are equivalent to the elasticities
$\varepsilon^{2}_{1},\varepsilon^{3}_{2}$, and $\varepsilon^{4}_{3}$
respectively.

MCA categories the logarithmic gains depending on their role and also
introduces additional logarithmic gains not present in BST although in
principle could be computed by BST. BST does have the additional sensitivity
measure that was mentioned above. Two of these logarithmic gains are relevant
here, they are the concentration control coefficients and concentration
response coefficients. The response coefficient describes how an external
factor such as $X_{0}$ influences the steady-state. The response coefficient
is equivalent to the logarithmic gain derived above using BST. In MCA
terminology the response coefficient of $X_{1}$ with respect to $X_{0}$ is
given by $R^{X_{1}}_{X_{0}}$:

$\displaystyle R^{X_{1}}_{X_{0}}=\frac{\partial\log X_{1}}{\partial\log
X_{0}}=L_{X_{1},X_{0}}=L_{ik}$ (10)

As noted before, un BST the logarithmic gains are also given the symbol
$L_{ik}$ which describes how variable $X_{i}$ is influenced by external factor
$k$.

In terms of the response coefficient A theorem in MCA states that a given
response coefficient is related to the product of an elasticity and a control
coefficient as follows:

$\displaystyle R^{X_{i}}_{X_{0}}=C^{X_{i}}_{v_{i}}\varepsilon^{v_{i}}_{X_{o}}$
(11)

where $C^{X_{i}}_{v_{i}}$ is the control coefficient of $X_{i}$ with respect
to reaction $v_{i}$. A concentration control coefficient describes how a
perturbation in a reaction (for example via a change in enzyme activity) can
lead to a change in a steady-state concentration. Thus $C^{X_{j}}_{v_{i}}$
describes how reaction $v_{i}$ influences species $X_{j}$.

In MCA the control coefficient for the first step, $v_{1}$, with respect to
$X_{1}$ can be shown to equal:

$\displaystyle C^{X_{1}}_{v_{1}}=\frac{1}{\varepsilon^{2}_{1}}$ (12)

This allows us to write the response coefficient as:

$\displaystyle
R^{X_{1}}_{X_{0}}=\frac{\varepsilon^{1}_{X_{o}}}{\varepsilon^{2}_{1}}$ (13)

which is again equivalent to the same equation from BST. The equivalence
between all the equations is summarized below:

$\displaystyle L_{X_{1},X_{0}}$
$\displaystyle=\frac{g_{10}}{h_{11}}=R^{X_{1}}_{X_{0}}=\frac{\varepsilon^{1}_{X_{o}}}{\varepsilon^{2}_{1}}$
(14) $\displaystyle L_{X_{2},X_{0}}$
$\displaystyle=\frac{g_{10}}{h_{22}}=R^{X_{2}}_{X_{0}}=\frac{\varepsilon^{1}_{X_{o}}}{\varepsilon^{3}_{2}}$
(15) $\displaystyle L_{X_{3},X_{0}}$
$\displaystyle=\frac{g_{10}}{h_{33}}=R^{X_{3}}_{X_{0}}=\frac{\varepsilon^{1}_{X_{o}}}{\varepsilon^{4}_{3}}$
(16)

### Including Product Inhibition

The previous example assumed that each reaction was unaffected by its product.
In reality, this is not necessarily the case (Shen et al., 2020). In this
example, we will include product inhibition. To keep things simple, let’s look
at a two-step pathway but it can be extended to any length pathway:

$X_{o}\stackrel{{\scriptstyle
v_{1}}}{{\longrightarrow}}X_{1}\stackrel{{\scriptstyle
v_{2}}}{{\longrightarrow}}$

If $X_{1}$ can product inhibit $v_{1}$, the power law for $v_{1}$ must be
modified to:

$v_{1}=\alpha_{1}X_{0}^{g_{10}}X_{1}^{h_{01}}$

The differential equation for $X_{1}$ then becomes;

$\frac{dX_{1}}{dt}=\alpha_{1}X_{0}^{g_{10}}X_{1}^{h_{01}}-\beta_{1}X_{1}^{h_{11}}$

Setting this equation to zero, taking logarithms on both sides and rearranging
to solve for $y_{1}$ as before yields:

$y_{1}=\frac{g_{10}y_{0}-b_{1}}{h_{11}-h_{01}}$

Compared to equation (5) we see there is an extra term in the denominator. As
before we can derive the logarithmic gain by differentiation to obtain:

$\frac{\partial\log X_{1}}{\partial\log X_{o}}=\frac{g_{10}}{h_{11}-h_{01}}$

If we use MCA, we obtain the same relationship:

$R^{X_{1}}_{X_{0}}=\frac{\varepsilon^{1}_{0}}{\varepsilon^{2}_{1}-\varepsilon^{1}_{1}}$

with the kinetic orders replaced by elasticities. Both BST and MCA produce the
same results.

### Including Negative Feedback

One final example worth looking at is when there is a negative feedback loop.
Here we will look at the results common to both as the BST community has done
a more extensive analysis of negative feedback. For comparison purposes, we
will use the same feedback model found in (Savageau, 1972). This model is a
four-step linear chain with a feedback loop from $X_{3}$ to the first step
(Figure 2).

Figure 2: Four-step pathway with negative feedback.

Following (Savageau, 1972) the power law for the inhibited step is given by:

$v_{1}=\alpha_{1}X_{0}^{g_{10}}X_{3}^{g_{13}}$

This is very similar to the way product inhibition was dealt with in the
previous section. Using the same approach as before, we write out the
differential equations, setting them to zero, and taking the log on both sides
yields:

$\displaystyle b_{1}$ $\displaystyle=g_{10}y_{0}-h_{11}y_{1}+g_{13}y_{3}$ (17)
$\displaystyle b_{2}$ $\displaystyle=h_{11}y_{1}-h_{22}y_{2}$ (18)
$\displaystyle b_{3}$ $\displaystyle=h_{22}y_{2}-h_{33}y_{3}$ (19)

Solving for the $y$ terms gives:

$\displaystyle y_{1}$ $\displaystyle=\frac{b_{1}\ h_{33}+b_{2}\ g_{13}+b_{3}\
g_{13}-g_{10}\ h_{33}\ y_{o}}{h_{11}(g_{13}-h_{33})}$ (20) $\displaystyle
y_{2}$ $\displaystyle=-\frac{h_{33}(b_{1}+b_{2}-g_{10}\ y_{o})+b_{3}\
g_{13}}{h_{22}\ (g_{13}-h_{33})}$ (21) $\displaystyle y_{3}$
$\displaystyle=\frac{b_{1}+b_{2}+b_{3}-g_{10}\ y_{o}}{g_{13}-h_{33}}$ (22)

In case you’re wondering, there is no $y_{o}$ solution because $y_{o}$ is the
input species which is fixed. To get the logarithmic gains we differentiate
these with respect to $y_{o}$ to yield (Savageau, 1972):

$\displaystyle\frac{\partial y_{1}}{\partial y_{0}}$
$\displaystyle=\frac{\partial\log X_{1}}{\partial\log
X_{o}}=\frac{h_{33}}{\left(h_{33}-g_{13}\right)}\frac{g_{10}}{h_{11}}$ (23)
$\displaystyle\frac{\partial y_{2}}{\partial y_{0}}$
$\displaystyle=\frac{\partial\log X_{2}}{\partial\log
X_{o}}=\frac{h_{33}}{\left(h_{33}-g_{13}\right)}\frac{g_{10}}{h_{22}}$ (24)
$\displaystyle\frac{\partial y_{3}}{\partial y_{0}}$
$\displaystyle=\frac{\partial\log X_{3}}{\partial\log
X_{o}}=\frac{h_{33}}{\left(h_{33}-g_{13}\right)}\frac{g_{10}}{h_{33}}$ (25)

The same analysis (see appendix) can be done using MCA to give:

$\displaystyle\begin{split}R^{X_{1}}_{X_{0}}&=\frac{\varepsilon^{4}_{3}}{(\varepsilon^{4}_{3}-\varepsilon^{1}_{3})}\frac{\varepsilon^{1}_{0}}{\varepsilon^{2}_{1}}\\\\[5.0pt]
R^{X_{2}}_{X_{0}}&=\frac{\varepsilon^{4}_{3}}{(\varepsilon^{4}_{3}-\varepsilon^{1}_{3})}\frac{\varepsilon^{1}_{0}}{\varepsilon^{3}_{2}}\\\\[5.0pt]
R^{X_{3}}_{X_{0}}&=\frac{\varepsilon^{4}_{3}}{(\varepsilon^{4}_{3}-\varepsilon^{1}_{3})}\frac{\varepsilon^{1}_{0}}{\varepsilon^{4}_{3}}\end{split}$
(26)

As we can see these equations are exactly the same as those derived via BST.

##  3 Discussion

This article discusses the development and comparison of two frameworks for
understanding biochemical networks: Biochemical Systems Theory (BST) and
Metabolic Control Analysis (MCA). BST, introduced by Savageau, focuses on
developing mathematical tools to analyze biochemical network behavior,
particularly emphasizing system stability and the use of power laws. MCA,
pioneered by Kacser and Burns, and Heinrich and Rapoport, approaches the same
problem but with a different emphasis, based on linearization of governing
equations and examining relationships between the various logarithmic gains.

There is a strong equivalence between BST and MCA despite their apparent
differences. This is shown through three simple examples of a linear chain of
biochemical reactions. Both frameworks offer similar is not identical insights
into how biochemical networks operate.

The key innovation of BST lies in its use of power laws to model biochemical
reactions, allowing for the straightforward formulation into systems of
differential equations. These equations facilitate the analysis of steady-
state concentrations and logarithmic gains, which represent how changes in a
given parameter affect variables such as concentrations in the system.

Comparing the derivations of logarithmic gains in BST and MCA, it becomes
evident that they yield equivalent results, albeit expressed differently. Both
frameworks provide a way to quantify the influence of external factors on
steady-state concentrations and highlight how local properties contribute to
system-wide behavior. Due to the way BST uses aggregation for more complex
networks, I feel working with more complex networks is better handled by MCA
because the entire procedure can be automated by computer (Christensen,
Hofmeyr and Rohwer, 2018) directly from an SBML model (Hucka et al., 2003).

## Acknowledgements

This work was partly supported by DOE award DESC0023091.

## Note to readers

I welcome critiques of this paper and can update the text accordingly.

## References

  * (1)
  * Alves et al. (2008) Alves, Rui, Ester Vilaprinyo, Benito Hernández-Bermejo and Albert Sorribas. 2008.  “Mathematical formalisms based on approximated kinetic representations for modeling genetic and metabolic pathways.” Biotechnology and Genetic Engineering Reviews 25(1):1–40.
  * Christensen, Hofmeyr and Rohwer (2018) Christensen, Carl D, Jan-Hendrik S Hofmeyr and Johann M Rohwer. 2018.  “PySCeSToolbox: a collection of metabolic pathway analysis tools.” Bioinformatics 34(1):124–125.
  * Heinrich and Rapoport (1974) Heinrich, Reinhart and Tom A Rapoport. 1974.  “A linear steady-state treatment of enzymatic chains: general properties, control and effector strength.” European journal of biochemistry 42(1):89–95.
  * Hucka et al. (2003) Hucka, Michael, Andrew Finney, Herbert M Sauro, Hamid Bolouri, John C Doyle, Hiroaki Kitano, Adam P Arkin, Benjamin J Bornstein, Dennis Bray, Athel Cornish-Bowden et al. 2003.  “The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models.” Bioinformatics 19(4):524–531.
  * Kacser (1973) Kacser, Henrik. 1973.  The control of flux. In Symp Soc Exp Biol.  Vol. 27 p. 65.
  * Kacser et al. (1995) Kacser, Henrik, James A Burns, H Kacser and DA Fell. 1995.  “The control of flux.” Biochemical Society Transactions 23(2):341–366.
  * Sauro (2018) Sauro, Herbert M. 2018.  Systems Biology: An Introduction to Metabolic Control Analysis.  Ambrosius Pubishing.
  * Savageau (1969) Savageau, Michael A. 1969.  “Biochemical systems analysis: II. The steady-state solutions for an n-pool system using a power-law approximation.” Journal of theoretical biology 25(3):370–379.
  * Savageau (1971) Savageau, Michael A. 1971.  “Concepts relating the behavior of biochemical systems to their underlying molecular properties.” Archives of Biochemistry and Biophysics 145(2):612–621.
  * Savageau (1972) Savageau, Michael A. 1972.  The behavior of intact biochemical control systems. In Current topics in cellular regulation.  Vol. 6 Elsevier pp. 63–130.
  * Savageau (1976) Savageau, Michael A. 1976.  Biochemical systems analysis: a study of function and design in molecular biology.  Addison-Wesley Pub.
  * Savageau (1988) Savageau, Michael A. 1988.  “Introduction to S-systems and the underlying power-law formalism.” Mathematical and Computer Modelling 11:546–551.
  * Shen et al. (2020) Shen, Lu, Martha Kohlhaas, Junichi Enoki, Roland Meier, Bernhard Schönenberger, Roland Wohlgemuth, Robert Kourist, Felix Niemeyer, David van Niekerk, Christopher Bräsen et al. 2020.  “A combined experimental and modelling approach for the Weimberg pathway optimisation.” Nature communications 11(1):1098.
  * Sorribas and Savageau (1989) Sorribas, Albert and Michael A Savageau. 1989.  “Strategies for representing metabolic pathways within biochemical systems theory: reversible pathways.” Mathematical Biosciences 94(2):239–269.
  * Voit (2013) Voit, Eberhard O. 2013.  “Biochemical systems theory: a review.” International Scholarly Research Notices 2013\.

## Appendix

One way to derive the control coefficients in MCA is to use implicit
differentiation of the steady-state equations. We start by writing the
differential equations in functional form:

$\displaystyle\frac{dX_{1}}{dt}$
$\displaystyle=v_{1}(X_{0})-v_{2}(X_{1}(X_{0}))$ (27)
$\displaystyle\frac{dX_{2}}{dt}$
$\displaystyle=v_{2}(X_{1}(X_{0}))-v_{3}(X_{2}(X_{0}))$ (28)
$\displaystyle\frac{dX_{3}}{dt}$
$\displaystyle=v_{3}(X_{3}(X_{0}))-v_{4}(X_{4}(X_{0}))$ (29)

For example, $v_{2}(X_{1}(X_{0}))$ means that the rate of reaction $v_{2}$ is
a function of substrate $X_{1}$. However, $X_{1}$ in turn is a function of the
fixed input substrate $(X_{0})$.

We set the equations to zero to indicate we are interested in the steady-
state:

$\displaystyle 0$ $\displaystyle=v_{1}(X_{0})-v_{2}(X_{1}(X_{0}))$ (30)
$\displaystyle 0$ $\displaystyle=v_{2}(X_{1}(X_{0}))-v_{3}(X_{2}(X_{0}))$ (31)
$\displaystyle 0$ $\displaystyle=v_{3}(X_{2}(X_{0}))-v_{4}(X_{3}(X_{0}))$ (32)

Since we are concerned with how $X_{0}$ influences the steady-state
concentrations, every steady-state species must be a function of $X_{0}$. To
obtain the logarithmic gains we differentiate each equation with respect to
$X_{0}$:

$\displaystyle 0$ $\displaystyle=\frac{\partial v_{1}}{\partial
X_{0}}-\frac{\partial v_{2}}{\partial X_{1}}\frac{dX_{1}}{dX_{0}}$ (33)
$\displaystyle 0$ $\displaystyle=\frac{\partial v_{2}}{\partial
X_{1}}\frac{dX_{1}}{dX_{0}}-\frac{\partial v_{3}}{\partial
X_{2}}\frac{dX_{2}}{dX_{0}}$ (34) $\displaystyle 0$
$\displaystyle=\frac{\partial v_{3}}{\partial
X_{2}}\frac{dX_{2}}{dX_{0}}-\frac{\partial v_{4}}{\partial
X_{3}}\frac{dX_{3}}{dX_{0}}$ (35)

These equations form a set of simultaneous equations that can be solved for
the three derivatives $dX_{1}/dX_{0},dX_{2}/dX_{0}$, and $dX_{3}/dX_{0}$. For
example, the first equation can be solved for $dX_{1}/dX_{0}$:

$\frac{dX_{1}}{dX_{0}}=\frac{\partial v_{1}}{\partial
X_{o}}\bigg{/}\frac{\partial v_{2}}{\partial X_{1}}$

The logarithmic gain can be obtained by scaling both sides, noting that at
steady-state $v_{1}=v_{2}$:

$\displaystyle\frac{dX_{1}}{dX_{0}}\frac{X_{0}}{X_{1}}=\frac{\partial
v_{1}}{\partial X_{0}}\frac{X_{0}}{v_{1}}\bigg{/}\frac{\partial
v_{2}}{\partial X_{1}}\frac{X_{1}}{v_{2}}$ (36)

The terms on the right-hand side are the elasticities, $\varepsilon^{1}_{0}$
and $\varepsilon^{2}_{1}$ respectively so that the equation can be written as:

$\displaystyle\frac{dX_{1}}{dX_{0}}\frac{X_{0}}{X_{1}}=\frac{\varepsilon^{1}_{0}}{\varepsilon^{2}_{1}}$
(37)

In MCA the left-hand term is called the response coefficient and designated
$R^{X_{1}}_{X_{0}}$. The remaining unknowns can be solved in the same way
leading to the equation shown in the main text (9).

Negative Feedback

Adding negative feedback adds a term to $v_{1}$ due to the negative feedback
via $X_{3}$. Compared to a system without negative feedback, only the first
differential equation is modified as shown below:

$\displaystyle\frac{dX_{1}}{dt}$
$\displaystyle=v_{1}(X_{0},X_{3}(X_{0}))-v_{2}(X_{1}(X_{0}))$ (38)
$\displaystyle\frac{dX_{2}}{dt}$
$\displaystyle=v_{2}(X_{1}(X_{0}))-v_{3}(X_{2}(X_{0}))$ (39)
$\displaystyle\frac{dX_{3}}{dt}$
$\displaystyle=v_{3}(X_{3}(X_{0}))-v_{4}(X_{4}(X_{0}))$ (40)

As before we set the equations to zero to find the steady-state:

$\displaystyle 0$
$\displaystyle=v_{1}(X_{0},X_{3}(X_{0}))-v_{2}(X_{1}(X_{0}))$ (41)
$\displaystyle 0$ $\displaystyle=v_{2}(X_{1}(X_{0}))-v_{3}(X_{2}(X_{0}))$ (42)
$\displaystyle 0$ $\displaystyle=v_{3}(X_{3}(X_{0}))-v_{4}(X_{4}(X_{0}))$ (43)

We implicitly differentiate with respect to $X_{0}$ to obtain:

$\displaystyle 0$ $\displaystyle=\frac{\partial v_{1}}{\partial
X_{0}}+\frac{\partial v_{1}}{\partial
X_{3}}\frac{dX_{3}}{dX_{0}}-\frac{\partial v_{2}}{\partial
X_{1}}\frac{dX_{1}}{dX_{0}}$ (44) $\displaystyle 0$
$\displaystyle=\frac{\partial v_{2}}{\partial
X_{1}}\frac{dX_{1}}{dX_{0}}-\frac{\partial v_{3}}{\partial
X_{2}}\frac{dX_{2}}{dX_{0}}$ (45) $\displaystyle 0$
$\displaystyle=\frac{\partial v_{3}}{\partial
X_{2}}\frac{dX_{2}}{dX_{0}}-\frac{\partial v_{4}}{\partial
X_{3}}\frac{dX_{3}}{dX_{0}}$ (46)

This forms a set of linear equations that can be solved for
$dX_{1}/dX_{0},dX_{2}/dX_{0}$, and $dX_{3}/dX_{0}$ to give:

$\displaystyle\frac{dX_{1}}{dX_{0}}$ $\displaystyle=\dfrac{\dfrac{\partial
v_{1}}{\partial X_{0}}\dfrac{\partial v_{4}}{\partial
X_{3}}}{\left(\dfrac{\partial v_{4}}{\partial X_{3}}-\dfrac{\partial
v_{1}}{\partial X_{3}}\right)\dfrac{\partial v_{2}}{\partial X_{1}}}$ (47)
$\displaystyle\frac{dX_{2}}{dX_{0}}$ $\displaystyle=\dfrac{\dfrac{\partial
v_{1}}{\partial X_{0}}\dfrac{\partial v_{4}}{\partial
X_{3}}}{\left(\dfrac{\partial v_{4}}{\partial X_{3}}-\dfrac{\partial
v_{1}}{\partial X_{3}}\right)\dfrac{\partial v_{3}}{\partial X_{2}}}$ (48)
$\displaystyle\frac{dX_{3}}{dX_{0}}$
$\displaystyle=\dfrac{1}{\left(\dfrac{\partial v_{4}}{\partial
X_{3}}-\dfrac{\partial v_{1}}{\partial X_{3}}\right)}$ (49)

These are unscaled expressions but it is straightforward to convert these to
logarithmic gains by multiplying both sides by $X_{0}$ and dividing both sides
by either $X_{1},X_{2}$ or $X_{3}$ depending on the equation. Because the
pathway is linear, at steady-state all rates are equal, thus
$v_{1}=v_{2}=v_{3}=v_{4}$ allowing us to scale the elasticities by $v$. The
result is the following scaled expressions as given in the main text (26):

$R^{X_{1}}_{X_{0}}=\frac{\varepsilon^{4}_{3}}{(\varepsilon^{4}_{3}-\varepsilon^{1}_{3})}\frac{\varepsilon^{1}_{0}}{\varepsilon^{2}_{1}}$

$R^{X_{2}}_{X_{0}}=\frac{\varepsilon^{4}_{3}}{(\varepsilon^{4}_{3}-\varepsilon^{1}_{3})}\frac{\varepsilon^{1}_{0}}{\varepsilon^{3}_{2}}$

$R^{X_{3}}_{X_{0}}=\frac{\varepsilon^{4}_{3}}{(\varepsilon^{4}_{3}-\varepsilon^{1}_{3})}\frac{\varepsilon^{1}_{0}}{\varepsilon^{4}_{3}}$


###### Abstract

An electric network whose dynamics is entirely equivalent to the Foucault
pendulum is presented. This circuit allows one to measure in a simulation of
the circuit, and eventually in a built circuit in the laboratory, the Hannay
geometric phase. This phase corresponds, via the equivalence discussed in the
text, exactly to the latitude effect on the rotation of the oscillation plane
of the pendulum. The circuit includes a gyrator, that introduces a non-
reciprocal character to the network.

MEASURING THE HANNAY GEOMETRIC PHASE

H. Fanchiotti1, C.A. García Canal1, M. Mayosky2, A. Veiga3, V. Vento4

1IFLP/CONICET and Departamento de Física
Universidad Nacional de La Plata, C.C.67, 1900, La Plata, Argentina
2LEICI, Departamento de Electrotecnia Facultad de Ingeniería
Universidad Nacional de La Plata, La Plata, Argentina and
Comisión de Investigaciones Científicas de la Provincia de Buenos Aires-
CICpBA, Argentina
3LEICI, Departamento de Electrotecnia Facultad de Ingeniería,Universidad
Nacional de La Plata, La Plata, Argentina; CONICET
4Departamento de Física Teórica-IFIC. Universidad de Valencia-CSIC.
E-46100, Burjassot (Valencia), Spain

The aim of this letter is to present a measurement of the geometric phase in a
classical system. This classical geometric phase, called Hannay angle [1] is
the analog to the well-known quantum geometric phase discovered by Berry [2].
We notice that in Ref.[3] (see also [4]) a particular circuit for obtaining
the Hannay angle was investigated based on time-dependent circuital elements
and a scheme for measuring this phase was proposed, wishing for an
experimental verification to be done. Another proposal was presented in
Ref.[5]. In this presentation we show an experimental electronic setup where
the mentioned phase can be measured in a simulation.

The standard example of a classical geometric phase is the latitude effect of
the Foucault pendulum. This pendulum allowed Leon Foucault in 1851 to
demonstrate the Earth’s rotation. The pendulum has been much studied and was
found that the oscillation plane of the pendulum has a rotation frequency that
depends on the latitude where it is located. In our opinion the Foucault
pendulum merits the representation in terms of an equivalent electronic
circuit that we present here because it provides an easy way of measuring the
classical geometric phase.

Let us start by remembering the standard equations that describe the Foucault
pendulum and that one can find in any textbook. Calling $x$ and $y$ the
projections of the bob of the pendulum in the plane $(x,y)$ perpendicular to
the vertical direction, the equations under the usual approximations are

$\ddot{x}+2\,\Omega_{T}\,\sin(\lambda)\,\dot{y}+\omega_{0}^{2}\,x=0,$ (1)
$\ddot{y}-2\,\Omega_{T}\,\sin(\lambda)\,\dot{x}+\omega_{0}^{2}\,y=0,$ (2)

where $\Omega_{T}$ is the Earth rotation frequency, $\lambda$ is the latitude
of the pendulum location and $\omega_{0}$ the oscillation frequency of the
pendulum.

G$C$$V_{C_{0}l}$$L$$V_{1}$$L$$V_{2}$$C$$V_{C_{0}r}$$t{=}0$$t{=}0$ Figure 1:
Resonant circuits coupled by a Gyrator

As is well known, an electric network includes a set of elements like
capacitors and inductances and connections among them. We are interested in
circuits where voltage and current depend only on time and we restrict our
analysis to passive networks. The network has ports: pairs of terminals that
allow to exchange energy with the surroundings and have a given voltage and
current. We call V the vector corresponding to the port voltage. A network not
connected to external energy sources is reciprocal if considering two
different terminals, say $\alpha$ and $\beta$, the excitation in $\alpha$
gives rise to a response in $\beta$ that is invariant under the permutation
$\alpha\Leftrightarrow\beta$.

In the area of electric networks one can find another classical system that is
governed by completely analog equations to those presented above for the
Foucault pendulum. It is the circuit network including two $LC$ oscillators
with the same frequency $\omega_{0}=1/\sqrt{LC}$, coupled by a so called
gyrator as presented in Fig.1. Due to the fact that any combination of $(L,C)$
elements gives rise to a reciprocal network, in order to have non-reciprocal
behavior the introduction of the gyrator, which is a passive element of two
ports[6] with conductance $G$, is needed.

The analysis of the time evolution of the electric circuit results in a system
of linear differential equations with constant coefficients that are [7]

$\ddot{V_{1}}+\gamma_{g}\,\dot{V_{2}}+\omega_{0}^{2}\,V_{1}=0,$ (3)
$\ddot{V_{2}}-\gamma_{g}\,\dot{V_{1}}+\omega_{0}^{2}\,V_{2}=0,$ (4)

where

$\gamma_{g}=\frac{G}{C}\,\,\,;\,\,\,{\omega_{0}}^{2}=\frac{1}{LC}.$

Consequently, the differential equations for both systems, Eqs.(1) and (2) and
Eqs. (3) and (4), have exactly the same form. In other words, the role of $x$
and $y$ is taken by $V_{1}$ and $V_{2}$ respectively, by establishing the
connection

$2\,\Omega_{T}\,\sin(\lambda)=\frac{G}{C}\,\,\,\,\,and\,\,\,\,\,\omega_{0}^{2}=\frac{1}{L\,C}.$

In the previous expressions, as we have stated before, $G$ represents the
conductivity of the gyrator while $C$ and $L$ are the capacity and the
inductance of each of the coupled oscillators. In order to obtain the
appropriate solutions of the differential equations they have to be
implemented with the corresponding initial conditions. We will return to this
point below.

Let us now go to discuss the appearence of geometric phase in our scheme. In
this respect, it is worth to remember that, as was mentioned, the motion of
the oscillation plane of the Foucault pendulum acquires a latitude effect in
the form of a phase shift [8]. This phase was shown [9] to be no more than an
example of the Hannay phase [1], the classical analog of the Berry phase. In
fact, whenever one performs the analysis of the pendulum motion as seen from
an inertial frame outside the Earth with the Earth’s center fixed, it is
possible to show [10] that the geometry of the phase is directly related to
the solid angle defined at the center of the Earth by the parallel of the
given latitude.

The analysis of the geometric phase in the case of the Foucault pendulum can
be presented as follows. With the projections $x$ and $y$ of the pendulum
movement we form the complex variable $z=x+iy$ that on the basis of equations
(1) and (2)) satisfies

$\ddot{z}+2\,i\,\Omega_{T}\sin(\lambda)\,\dot{z}+a^{2}\,z=0,$ (5)

whose general solution reads

$z(t)=z_{0}(t)e^{-i\Omega_{T}\,\sin(\lambda)\,t},$ (6)

with

$z_{0}(t)=A_{1}e^{ia\,t}+A_{2}e^{-ia\,t},$

and $A_{1}$ and $A_{2}$ are constants fixed by the initial conditions.

Then, at $t=0$ one has

$z(t)=z_{0}(0)\,\,\,;\,\,\,z_{0}(0)=A_{1}+A_{2},$

and we can evaluate the solution for a time $t=T$ the period of rotation of
the Earth, namely

$T=\frac{2\pi}{\Omega_{T}},$

to obtain

$z(T)=z_{0}(T)e^{-i2\pi\sin(\lambda)}\,\,\,;\,\,\,z_{0}(T)=A_{1}e^{i2\pi\omega_{0}/\Omega_{T}}+A_{2}e^{-i2\pi\omega_{0}/\Omega_{T}}.$

The pendulum frequency is so large when compared with $\Omega_{T}$, that one
can safely consider

$\omega_{0}/\Omega_{T}=n.$ (7)

Consequently, the solution has the expression

$z(T)=z_{0}(0)e^{-i2\pi\sin(\lambda)},$ (8)

that clearly shows that the solution $z$ returns to the initial state with an
extra phase given by

$\eta=2\,\pi\,\sin(\lambda).$ (9)

This is the Hannay phase associated to the latitude effect of this classical
system. This effect is described by the angle $\Delta\phi$ that the
oscillation plane of the pendulum has to travel after a complete rotation of
the Earth to complete a full $2\pi$ rotation in a given latitude which is
given by

$\triangle\phi=2\,\pi-\eta=2\,\pi(1-cos(\lambda)).$ (10)

Let us now proceed with the analysis of the circuit of Fig.1.

### Results of the simulation

As we stated, the electric equivalence comprises two identical resonant
circuits connected by a gyrator, as shown in Fig.1. The behavior of this
circuit can be simulated using standard packages such as LT-spice. Although it
is not included as a component in the LT-spice standard libraries, an ideal
gyrator can be implemented with two voltage dependent current sources in a
straightforward way. The plot of $V_{1}$ versus $V_{2}$ resembles exactly the
time evolution of the $x$ and $y$ projections of the pendulum bob. Initial
conditions for simulation are imposed as initial charges in each capacitor.
This guarantees that this initial conditions represent the Foucault pendulum
when it starts with zero initial velocity.

The geometric phase equivalent is now

$\eta=\frac{\pi}{\Omega_{T}C}G,$ (11)

so that different latitudes can be simulated just by changing the value of $G$
multiplying it times a value in between $0$ and $1$. In Fig.2, we show the
results of the simulations that correspond in the case of the pendulum to
latitudes of $90^{\circ}$ (pole) and $60^{\circ}$. Notice that for clarity we
decided to show only half of the $2\,\pi$ rotation (half of a day). For these
simulations, the circuit values are the following:

$\displaystyle L$ $\displaystyle=$ $\displaystyle 10mHy$ $\displaystyle C$
$\displaystyle=$ $\displaystyle 0.1\mu F$ $\displaystyle G$ $\displaystyle=$
$\displaystyle 52.354\mu S*\sin(latitude)$

Initial conditions for both cases are $V_{C_{0}l}=1V$ and $V_{C_{0}r}=0V$.
Initial currents on inductors are both null. For these circuit values a 12
hour pendulum movement corresponds to 12 milliseconds of simulation time.

Figure 2: Simulation results. (a) latitude $90^{\circ}$. (b) latitude
$60^{\circ}$

It is worth mentioning that in an experimental circuit implementation
unavoidable component losses arise. In such a case, the dissipation due to the
resistances involved should certainly be taken into account. This implies the
modification of the equations (3) and (4), and both the simulation and the
experimental implementation get more involved. For pedagogical purposes, we
postpone the treatment of these complexities to a forthcoming publication.

In conclusion, the measurement of the geometric phase of a classical system,
the Hannay angle, is easily realized in a simulation of an electric network
including a gyrator. This determination has direct correspondence with the so-
called latitude effect in the rotation of the oscillation plane of the
Foucault pendulum. To summarize, besides having described a classical system
where the geometric phase can be measured, we have presented a complete
electric analog of the Foucault pendulum and have made a simulation that
allows to measure geometric angles.

## Acknowledgement

This work was supported in part by the MICINN and UE Feder under contract
FPA2016-77177-C2-1-P.

## References

  * [1] J. H. Hannay, ”Angle variable holonomy in adiabatic excursion of an integrable Hamiltonian”, J. Phys. A 18, 221 (1985).
  * [2] M. V. Berry,”Quantal Phase Factors Accompanying Adiabatic Changes”, Proc. Roy. Soc. Lond. A 392 (1802), 45 (1984).
  * [3] D. Xu, ”Hannay angle in an LCR circuit with time-dependent inductance, capacity and resistance”, J. Phys. A:Math Gen. 35 L455 (2002).
  * [4] J.H. Hannay, ”Comment on ’Hannay angle in an LCR circuit with time-dependent inductance, capacity and resistance”, J. Phys. A:Math Gen. 35 9699 (2002).
  * [5] Sharba Bhattacharjee et al,”Study of geometric phase using classical coupled oscillators”, Eur. J. Phys. 39, 035404 (2018)
  * [6] N. Balabanian, T.A. Bickart, ”Linear network theory: analysis, properties, design and synthesis”, Willey (1969); H. Carlin, A. Giordano, ”Network Theory: An Introduction to Reciprocal and Nonreciprocal Circuits”, Prentice Hall (1964)
  * [7] M. Caruso, H. Fanchiotti, C.A. Garcia Canal, M. Mayosky, A. Veiga, ”The quantum CP-violating kaon system reproduced in the electronic laboratory”, Proc.Roy.Soc.Lond.A 472, 2195, 20160615 (2016)
  * [8] W. B. Somerville, ”The Description of Foucault’s Pendulum”, — Q. Jl. Astr. Soc. 13, 40 (1972).
  * [9] J. B. Hart, R. E. Miller and R. L. Miller, ”A simple geometric model for visualizing the motion of a Foucault pendulum”, Am. J. Phys. 55, 67 (1987); A. Khein and D. F. Nelson, ”Hannay angle study of the Foucault pendulum in action-angle variables”, Am. J. Phys. 61, 170 (1993).
  * [10] J. von Bergmann, ”Foucault pendulum through basic geometry”, Am. J. Phys. 75, 888 (2007)


# PerfVis: Pervasive Visualization in Immersive Augmented Reality for
Performance Awareness

Leonel Merino  VISUS, University of Stuttgart ,  Mario Hess  SCG, University
of Bern ,  Alexandre Bergel  ISCLab, DCC, University of Chile ,  Oscar
Nierstrasz  SCG, University of Bern and  Daniel Weiskopf  VISUS, University of
Stuttgart

(2019)

###### Abstract.

Developers are usually unaware of the impact of code changes to the
performance of software systems. Although developers can analyze the
performance of a system by executing, for instance, a performance test to
compare the performance of two consecutive versions of the system, changing
from a programming task to a testing task would disrupt the development flow.
In this paper, we propose the use of a city visualization that dynamically
provides developers with a pervasive view of the continuous performance of a
system. We use an immersive augmented reality device (Microsoft HoloLens) to
display our visualization and extend the integrated development environment on
a computer screen to use the physical space. We report on technical details of
the design and implementation of our visualization tool, and discuss early
feedback that we collected of its usability. Our investigation explores a new
visual metaphor to support the exploration and analysis of possibly very large
and multidimensional performance data. Our initial result indicates that the
city metaphor can be adequate to analyze dynamic performance data on a large
and non-trivial software system.

††copyright: rightsretained††journalyear: 2019††copyright:
acmlicensed††conference: Tenth ACM/SPEC International Conference on
Performance Engineering Companion; April 7–11, 2019; Mumbai, India††booktitle:
Tenth ACM/SPEC International Conference on Performance Engineering Companion
(ICPE ’19 Companion), April 7–11, 2019, Mumbai, India††price: 15.00††doi:
10.1145/3302541.3313104††isbn: 978-1-4503-6286-3/19/04††ccs: General and
reference Performance††ccs: Human-centered computing Visualization design and
evaluation methods††ccs: Human-centered computing Mobile devices††ccs:
Software and its engineering Software maintenance tools

##  1\. Introduction

Among the many questions that arise during software development, programmers
often formulate questions about the performance of software systems (LaToza
and Myers, 2010; Alcocer et al., 2013). They ask, for instance, “what is the
software doing when performance issues arise?” (Reiss and Renieris, 2005), and
“where is most of the time being spent?” (Shah and Guyer, 2016). Since
developers introduce multiple changes to the source code during the
implementation of a software system, being aware of the impact of a code
change task to the performance of a software system is an important concern
for developers (LaToza and Myers, 2010). Tools such as profilers and
performance tests allow developers to analyze the performance of only a single
version of a system at a time, and force developers to change to a different
task (interrupting their flow). Our intuition is that developers would benefit
from analyzing the live performance of an evolving system as its source code
changes. We expect that such an analysis would help them identify changes that
severely decrease the performance of a system.

Figure 1. Pervasive visualization displayed in immersive augmented reality: ⓐ
a developer wears a Microsoft HoloLens device ⓑ for performance visualizations
using the city metaphor and a scatter plot.

We conjecture that a pervasive tool, (_i.e._ , a tool that is omnipresent at
all times during the implementation of a software system) can make developers
aware of important changes in the performance of a system without disrupting
the implementation flow. Using a visualization approach, developers could
obtain such a pervasive view of the live performance of a system, however, we
observe that developers can be reluctant to sacrifice space in their
integrated development environment (IDE). Consequently, we formulate the
following research question:  _RQ: How can visualization support developers in
the analysis of the impact of source code changes to the performance of a
system?_

This paper presents an initial step toward answering our general research
question. To increase developer awareness of the impact of changes to the
source code in the performance of software systems we present _PerfVis_ , a
pervasive visualization tool displaying in immersive augmented reality.
PerfVis provides an omnipresent overview of the system that takes advantage of
visualization techniques to make developers aware of the impact of changes
introduced to the source code in the live performance of a software system.

The main contribution of the paper is threefold:  (i) a discussion of design
choices,  (ii) lessons learned from building the tool, and  (iii) insights
from observing its use.  We also contribute to the reproducibility of our
research by making the implementation of PerfVis publicly
available111http://scg.unibe.ch/research/perfvis.

##  2\. Related work

A number of 3D visualization tools have been proposed to support software
development concerns. _SeeIt 3D_ (Sharif et al., 2013) uses a 3D scatter plot
visualization to represent software metrics collected from Java source code to
support developers on software comprehension tasks. _SynchroVis_ (Waller et
al., 2013) employs a visualization based on the 3D city metaphor that displays
the structure and metrics of software systems, as well as program traces and
synchronization aspects (_e.g._ , semaphores/wait) to support the analysis of
the behavior of concurrent software systems. These tools display the 3D
visualization on a standard computer screen and focus on performance data
collected using static analysis. In contrast, we propose the visualization to
be displayed in immersive augmented reality and focus on live performance data
collected using dynamic analysis.

There are some visualization tools that use virtual reality to support
software comprehension tasks. _Code Park_ (Khaloo et al., 2017) proposes an
immersive visualization of software metrics and source code; an evaluation
showed that the visualization displayed in virtual reality excels at usability
and significantly helps in code understanding. _ExplorViz_ (Fittkau et al.,
2017) introduces a visualization of software landscapes to support software
comprehension tasks. _CityVR_ (Merino et al., 2017c) gamifies software
implementations using a 3D city visualization displayed in virtual reality to
boost developer engagement. In contrast, we propose the use of immersive
augmented reality to improve the awareness of the impact of changes on the
source code to the performance of software systems. Moreover, some of these
tools focus on various static software metrics, as opposed to our tool that
includes the visualization of dynamic live performance data.

Several existing performance visualizations focus on the analysis of static
performance data. _JOVE_ (Reiss and Renieris, 2005) is a visualization tool
for monitoring the performance of Java programs. Using the visualization,
developers can obtain details on demand to identify routines and threads in
which the application spends much time. Another tool (Blochinger et al., 2005)
includes the visualization of the structure of parallel software systems. The
visualization uses an execution graph to simplify the analysis of the complex
run-time behavior. Moreta and Telea (Moreta and Telea, 2007) present a
visualization of the behavior of the memory allocator in C programs to
optimize functionality, decrease fragmentation, and improve response time.
Another tool (Ogami et al., 2017) facilitates the visualization of software
performance in real-time, using the city metaphor to show the structure and
performance of a program. The buildings in the city represent the classes in
the system, and the heights of the buildings represent the number of times the
methods of a class are called during execution. All these tools support the
analysis of various aspects of software performance through visualizations
displayed on a standard computer screen. In contrast, we rely on immersive
augmented reality to display our visualization.

A few software visualizations chose augmented reality as a display medium. One
study (Merino et al., 2018) evaluated a city visualization displayed in
immersive augmented reality to support software comprehension tasks. Another
study (Souza et al., 2012) used the city metaphor in augmented reality to
support the analysis of software evolution. In contrast, we propose a city-
based visualization using immersive augmented reality to support software
performance tasks, which is an innovative and not-yet explored approach.

##  3\. PerfVis Overview

PerfVis is a tool to visualize software performance through immersive
augmented reality. We characterize the context of PerfVis using a taxonomy
proposed by Maletic _et al._ (Maletic et al., 2002) and Merino _et al._
(Merino et al., 2017b).  _PerfVis_ provides visualization in _immersive
augmented reality_ to support the _programmer_ audience in system performance
tasks. The visualization shows _static data_ of the structure of the system
and _dynamic data_ of the live performance of the system.

Figure 2. A detailed view of the ⓐ dynamic city visualization and ⓑ the
complementary scatter plot that maintains the history of the overall
performance of the system for a configurable interval of time.

###  3.1. Design

Medium: Immersive Augmented Reality. Most visualizations that support
performance analysis tasks are displayed on a computer screen. They require
developers to either use an application outside their development environment,
or use a plug-in in their environment that sacrifices some valuable screen
space. Amongst the few visualizations that use a different medium most are
displayed in immersive virtual reality (Merino et al., 2017a). These
visualizations require that developers completely change their development
environment, and isolate themselves wearing a headset. As opposed to these
approaches, we argue that a medium that disrupts as little as possible the
development flow can significantly increase the usability of visualizations.
The medium used to display a software visualization can impact its
effectiveness and efficiency (Merino, 2018). Consequently, we chose immersive
augmented reality to display PerfVis, because  (i) developers can use the
visualization as a complement to their usual development environment without
much disruption (_i.e._ , developers are required to wear a headset, but not
to modify the rest of their environment), and  (ii) visualizations in
immersive augmented reality are inherently displayed in 3D, which provides an
extra dimension to encode performance data of a software system.

Apparatus: HoloLens. We choose Microsoft HoloLens to render our visualization
in immersive augmented reality. HoloLens is an untethered head-mounted
stereoscopic display, featuring two stereoscopic displays with 1268$\times$720
pixel resolution, 60 Hz content refresh rate, and a 30° H and 17.5° V field of
view. Interaction is possible through head tracking, gesture input, and voice
support.

Technique: City Visualization and Scatter plot. In a previous study (Merino et
al., 2017b), we observed that most visualizations to support developers in
monitoring performance used pixel-based techniques. These visualizations
usually attempt to provide an overview of a significant amount of data to
promote the detection of visual patterns that may signal anomalies in the
performance of software systems. Instead, we adopt lightweight visualization
techniques that involve less cognitive effort to make developers aware of the
performance of a system. Our visualization includes two views:  (i) a city
visualization that provides a pervasive overview of the structure and
performance and  (ii) a scatter plot that displays the history of the
performance of a software system in the past.

_City visualization._ We opted to represent a metric of the performance of the
software system using the color of buildings. The visualization is always on
and indicates the system performance collected during a short time frame. In
consequence, a glimpse at the city can make developers aware of the current
overall status of the performance, and a detail of the classes that are
involved. Figure 2 \smalla⃝ shows a visualized software city. Each building
represents a class in the software system. Buildings are grouped into
districts according to their packages. Each building can be configured to
represent software metrics using three properties:  (i) height,  (ii)
footprint (_e.g._ , squared base), and  (iii) color.  We consequently
configured the city in the following manner: each building’s  (i) height
encodes the number of methods of the represented class,  (ii) footprint
(_i.e._ , squared base) encodes the number of attributes of the represented
class, and  (iii) color encodes the number of times that methods of the class
are called during execution.

_Scatter plot._ The scatter plot helps developers analyze changes in the
dynamic city visualization that might be overlooked. Figure 2 \smallb⃝ shows
the scatter plot to the right of the city visualization. The $X$-axis encodes
time, and the $Y$-axis encodes the classes in the system. Each mark represents
a class that was involved in the execution of the system at a point in time.
The color of a mark encodes the number of times methods of a class were
called. The more intense the red, the more often a class was called. Each row
in the scatter plot represents the evolution of the number of times the
methods of a class are called along time. Each column is a snapshot of the
performance of the whole system at a point in time.

Interaction: Selection and Navigation. Developers can interact with PerfVis
through:  (i) _Selection_. Through head movements, users can point and select
elements in the visualization. When the pointer hovers over a building in the
city visualization or a mark in the scatter plot, the name of the
corresponding class is displayed. Users select buildings by performing an
_airtap_ gesture.222https://docs.microsoft.com/en-us/hololens/ Once a building
is selected, that building and the mark that represent the same class are
simultaneously highlighted with a yellow background. Similarly, when selecting
a mark in the scatter plot, the building and mark that represent the same
class are highlighted. Using this feature, developers relate the information
they get from one visualization to the elements in the complementary
visualization.  (ii) _Navigation_. To gain an overview of the whole system,
the user can move around and inspect the city from different angles. This way
the user can compare metrics between different buildings.

###  3.2. Implementation

The landscape of the city visualization is generated using
CodeCity333http://smalltalkhub.com/#!/~RichardWettel/CodeCity for Moose
5.444http://www.moosetechnology.org/ The color animation of the city and the
scatter plot are implemented in Unity555Version 2017.3.0b8 (64-bit) using the
MixedRealityToolkit-Unity666https://github.com/Microsoft/MixedRealityToolkit-
Unity library. The performance data are collected using the
Spy2777https://github.com/ObjectProfile/Spy2 profiler for Pharo Smalltalk, and
streamed to the HoloLens through the network.

##  4\. Discussion

We revisit our research question: _How can visualization support developers in
the analysis of the impact of source code changes to the performance of a
system?_ We have presented a design and publicly available implementation of a
visualization tool that uses immersive augmented reality to support developers
in the analysis of software performance. We now discuss lessons learned during
the development process, and early feedback on the usability of PerfVis.

###  4.1. Lessons Learned

Scale. PerfVis aims to give developers an overview of the live performance of
a software system with the least disruption possible to the integrated
development environment. We tested PerfVis visualizations at various scales.
We observed that the visualization had to be scaled to be large enough to give
an overview of the complete system, but also small enough to be used as a
complement to a computer screen.

Text. We learned that the position, rotation, and scale of text can greatly
impact usability when navigation is done by head movements. In an early
prototype, we chose to display the names of classes on top of their
corresponding building, however, we found that depending on the position of
the user and the size of the building, the text on top of a building might be
difficult to read, and also that text on top of a building might be occluded
by neighboring buildings. We addressed these issues by positioning the text in
relation to the position of the camera. Specifically, we placed the text at
the top left corner of the user’s field of view. Benefits of this approach are
that  (i) text is never occluded,  (ii) text always faces the user, and  (iii)
no movement of the head and only little movement of the eyes is required to
read the text.

History. A live visualization has the benefit that it can provide developers
online feedback of the performance of the system, however, changes in the
visual properties of the visualization can take a very short time. Then the
city visualization can be good at making developers aware of an important
change in the performance but might not be good for analyzing details of the
context of the change since it provides developers too little time to react.
To help developers get details on demand we included a companion lightweight
visualization (_i.e._ , scatter plot). We observed some usability issues. For
instance,  (i) marks in the scatter plot can be too small and difficult to
select, and  (ii) the size of the buffer presented in the scatter plot has to
ensure a good performance (too many marks would slow down the visualization
tool).  To address the former issue we included a type of selection that
involves only head movement (marks are selected when hovering over them for
one second). We also addressed the latter issue by making the definition of
the size of the buffer configurable.

###  4.2. Early Feedback on Usability

We demonstrated the tool to a group of Computer Science students to discuss
our design choices and obtain impressions. Their feedback allowed us to
iterate several times on the design of the various components of the
visualization tool. Certainly, the positive impressions of developers who used
our visualization tool do not represent a solid proof of its benefits,
however, they represent a basis to design a suitable evaluation (Merino et
al., 2018).

Visualization Techniques. Developers who used PerfVis noticed that the city
visualization and the scatter plot complement each other. The former, being a
space-filling technique, can be scaled down to a small size, and so provide
developers an omnipresent overview of the performance of the whole system; and
once the attention of a developer is caught, the latter can be used to provide
details on demand of the history of the performance of a system. Some
developers who used early prototypes of PerfVis observed that the city
visualization was too big, and so only part of it was visible through the
narrow field of view of the immersive augmented reality device that we used
(HoloLens). We then scaled the visualization to be visible even when the
developer is looking at the computer screen.

Interactions. Interactions such as navigation and selection represent key
features that promote usability. In particular, developers who used PerfVis
found it difficult to stabilize their head when pointing at small elements in
the visualization to select them (through airtap). Some of the developers also
said that the highlighting feature of selected elements in the city and
scatter plot visualizations did not have a good contrast, and it was hard to
see. Several developers suggested including more interactions to the scatter
plot visualization. Specifically, they said that features to pause and rewind
would help them in the analysis of the performance of the system in a past
time.

##  5\. Conclusion and Future work

We introduced _PerfVis_ , a visualization tool displayed in immersive
augmented reality that supports developers in the analysis of software
performance. We elaborated on design choices and discussed implementation
details, lessons learned, and early impressions of developers who used
prototypes of our visualization tool. In the future, we plan to conduct a
thorough evaluation of PerfVis to identify strengths and improvement
opportunities.

## Acknowledgments

Merino and Weiskopf acknowledge funding by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) – Projektnummer 251654672 – TRR 161.
Nierstrasz thanks the Swiss National Science Foundation for its financial
support of “Agile Software Analysis” (project 162352). Bergel thanks LAM
Research for its financial support.

## References

  * (1)
  * Alcocer et al. (2013) Juan-Pablo Sandoval Alcocer, Alexandre Bergel, Stéphane Ducasse, and Marcus Denker. 2013.  Performance evolution blueprint: Understanding the impact of software evolution on performance. In Proc. of VISSOFT. IEEE, 1–9.
  * Blochinger et al. (2005) Wolfgang Blochinger, Michael Kaufmann, and Martin Siebenhaller. 2005.  Visualizing structural properties of irregular parallel computations. In Proc. of SOFTVIS. ACM, 125–134.
  * Fittkau et al. (2017) Florian Fittkau, Alexander Krause, and Wilhelm Hasselbring. 2017\.  Software landscape and application visualization for system comprehension with ExplorViz.  Information and Software Technology 87 (2017), 259–277.
  * Khaloo et al. (2017) Pooya Khaloo, Mehran Maghoumi, Eugene Taranta, David Bettner, and Joseph Laviola. 2017.  Code Park: A new 3D code visualization tool. In Proc. of VISSOFT. IEEE, 43–53.
  * LaToza and Myers (2010) Thomas D. LaToza and Brad A. Myers. 2010.  Hard-to-answer questions about code. In Proc of. PLATEAU. ACM, Article 8, 6 pages.
  * Maletic et al. (2002) Jonathan I. Maletic, Andrian Marcus, and Michael Collard. 2002\.  A task oriented view of software visualization. In Proc. of VISSOFT. IEEE, 32–40.
  * Merino (2018) Leonel Merino. 2018\.  The Medium of Visualization for Software Comprehension.  PhD thesis. University of Bern.
  * Merino et al. (2018) Leonel Merino, Alexandre Bergel, and Oscar Nierstrasz. 2018\.  Overcoming issues of 3D software visualization through immersive augmented reality. In Proc. of VISSOFT. IEEE, 54–64.
  * Merino et al. (2017a) Leonel Merino, Johannes Fuchs, Michael Blumenschein, Craig Anslow, Mohammad Ghafari, Oscar Nierstrasz, Michael Behrisch, and Daniel Keim. 2017a.  On the impact of the medium in the effectiveness of 3D software visualization. In Proc. of VISSOFT. IEEE, 11–21.
  * Merino et al. (2017c) Leonel Merino, Mohammad Ghafari, Craig Anslow, and Oscar Nierstrasz. 2017c.  CityVR: Gameful software visualization. In Proc. of ICSME (TD Track). IEEE, 633–637.
  * Merino et al. (2018) Leonel Merino, Mohammad Ghafari, Craig Anslow, and Oscar Nierstrasz. 2018.  A systematic literature review of software visualization evaluation.  Journal of Systems and Software 144 (2018), 165–180.
  * Merino et al. (2017b) Leonel Merino, Mohammad Ghafari, and Oscar Nierstrasz. 2017b.  Towards actionable visualization for software developers.  Journal of Software: Evolution and Process 30, 2 (2017), e1923.
  * Moreta and Telea (2007) Sergio Moreta and Alexandru Telea. 2007.  Visualizing dynamic memory allocations. In Proc. of VISSOFT. IEEE, 31–38.
  * Ogami et al. (2017) K. Ogami, R. G. Kula, H. Hata, T. Ishio, and K. Matsumoto. 2017.  Using high-rising cities to visualize performance in real-time. In Proc. of VISSOFT. IEEE, 33–42.
  * Reiss and Renieris (2005) Steven P. Reiss and Manos Renieris. 2005.  JOVE: Java as it happens. In Proc. of SOFTVIS. ACM, 115–124.
  * Shah and Guyer (2016) Michael D. Shah and Samuel Z. Guyer. 2016.  An interactive microarray call-graph visualization. In Proc. of VISSOFT. IEEE, 86–90.
  * Sharif et al. (2013) Bonita Sharif, Grace Jetty, Jairo Aponte, and Esteban Parra. 2013.  An empirical study assessing the effect of SeeIT 3D on comprehension. In Proc. of VISSOFT. IEEE, 1–10.
  * Souza et al. (2012) Rodrigo Souza, Bruno Silva, Thiago Mendes, and Manoel Mendonça. 2012.  SkyscrapAR: An augmented reality visualization for software evolution. In Proc. of WBVS. SBC, 17–24.
  * Waller et al. (2013) Jan Waller, Christian Wulf, Florian Fittkau, Philipp Dohring, and Wilhelm Hasselbring. 2013\.  Synchrovis: 3D visualization of monitoring traces in the city metaphor for analyzing concurrency. In Proc. of VISSOFT. IEEE, 1–4.