# BBP-type formulas — an elementary approach

SIMON KRISTENSEN  S. Kristensen, Department of Mathematics, Aarhus University,
Ny Munkegade 118, DK-8000 Aarhus C, Denmark<EMAIL_ADDRESS>and  OSKAR
MATHIASEN  O. Mathiasen, Department of Mathematics, Aarhus University, Ny
Munkegade 118, DK-8000 Aarhus C, Denmark<EMAIL_ADDRESS>
###### Abstract.

We provide a simple way of searching for formulas of the
Bailey–Borwein–Plouffe type together with an algorithm and an implementation
in sage. Aside from rediscovering some already known formulas, the method has
been used in the discovery of a new BBP-type formula for $\sqrt{3}\pi$. In
addition, the implementation is very flexible and allows us to look for BBP-
type formulas to irrational bases but with integer coefficients. As an example
of this, searching in various Pisot bases, we have discovered a formula for
$\pi$ in base $1+\sqrt{3}$, along with additional formulas.

###### Key words and phrases:

BBP-type formulas; null formulas; formulas for $\pi$

###### 1991 Mathematics Subject Classification:

Primary 11K16; Secondary 11Y60

##  1\. Introduction

In 1997, D. Bailey, P. B. Borwein and S. Plouffe [4] discovered the formula

(1)
$\pi=\sum_{k=0}^{\infty}\frac{1}{16^{k}}\left(\frac{4}{8k+1}+\frac{-2}{8k+4}+\frac{-1}{8k+5}+\frac{-1}{8k+6}\right)$

The novelty of the formula is that it allows one to calculate – using a spigot
algorithm – the $n$’th digit in the hexadecimal expansion of $\pi$ without
knowing the preceeding digits. Since then, many other formulas of similar type
have been discovered. We refer the reader to D. Bayley’s on-line compendium
[3] for a list of some of the known formulas.

Before proceeding, we will need some notation.

######  Definition 1.

A BBP-type formula is a series of the form

$\operatorname{BBP}(d,b,n,A)=\sum_{k=0}^{\infty}\frac{1}{b^{k}}\sum_{j=1}^{n}\frac{a_{j}}{(kn+j)^{d}},$

where $d,b,n\in\mathbb{N}$ are called the degree, base and number
respectively, and $A=(a_{1},\dots,a_{n})\in\mathbb{R}^{n}$ is a vector.

Often, one requires the $a_{j}$’s to be integers, and clearly this is the more
interesting case. However, we will not be making this restriction in
calculations, although we will present only results with integer coefficients.
In this way, (1) states that

$\pi=\operatorname{BBP}(1,16,8,(4,0,0,-2,-1,-1,0,0)).$

A unified approach to searching for BBP-type formulas was however sadly
lacking until D. Barsky, V. Muñoz and R. Pérez-Marco [5] provided a framework
for searching for such formulas. Their idea is to set up a family of iterated
integrals, depending on a single complex variable, $s$ say. These integrals
are subsequently expanded into rapidly converging series. On substituting
particular values of $s$ into the series and considering real and imaginary
parts separately, one obtains BBP-type formulas for a variety of numbers.

A main feature of their work is that they find representations of zero as a
BBP-type formula. These so-called null formulas have been among the more
mysterious ones in literature, and have given rise to relations between BBP-
type formulas.

In this paper, we give a less general form of the main theoretical tool of
Barsky, Muñoz and Pérez-Marco [5], which as it turns out admits a simple and
intuitive geometric interpretation. In addition, it is rather simple to
implement an algorithm based on the main result, which is flexible enough to
deal with a number of interesting cases.

##  2\. Main theoretical results

We present here our main theoretical result, which in its simplicity provides
several BBP-type formulas.

######  Theorem 1.

Let $\theta=\frac{a}{b}2\pi$ with $\frac{a}{b}\in\mathbb{Q}$ and $2|b$ and let
$r\in(0,1)$. Then,

$\arg(1+r\cdot
e^{i\theta})=\operatorname{BBP}\left(1,\tfrac{1}{r^{b}},b,A\right),$

where $A=(a_{1},\dots,a_{b})$ with $a_{j}=r^{j}(-1)^{j+1}\sin(j\theta)$.

###### Proof.

Our starting point is the Taylor series for $\log(1+z)$ centred at $z=0$,

$\log(1+z)=\sum_{n=1}^{\infty}(-1)^{n-1}\frac{z^{n}}{n}.$

Substituting $z=re^{i\theta}$, we get,

$\log(1+re^{i\theta})=\sum_{n=1}^{\infty}(-1)^{n-1}\frac{r^{n}e^{in\theta}}{n}.$

Recalling that the argument of a complex number is nothing but the imaginary
part of the logarithm of the number,

$\arg(1+re^{i\theta})=\sum_{n=1}^{\infty}(-1)^{n-1}\frac{r^{n}\sin(n\theta)}{n}.$

Now, $\sin(n\theta)$ is clearly $b$-periodic, so writing each $n=bk+j$ with
$j=1,2,\dots,b$, we get since $2|b$,

$\arg(1+re^{i\theta})=\sum_{j=0}^{\infty}\left(\frac{1}{r^{b}}\right)^{-k}\sum_{j=1}^{b}\frac{r^{j}(-1)^{j+1}\sin(j\theta)}{bk+j},$

as required. ∎

Note that we could just as well have taken the real part in the proof to
obtain a BBP-type formula for $\log|1+r\cdot
e^{i\theta}|=\frac{1}{2}\log(1+r^{2}+2r\cos\theta)$. In the example given
below, this would have resulted in a BBP-type formula for $\log 2$. The same
procedure gives a formula for the logarithm of a certain real number for each
of our obtained formulas below. We have left it for the interested reader to
deduce such formulas from the new ones obtained in this paper.

As an example of our main theorem, let us derive a BBP-type formula for
$\pi/4$, originally found by Ferguson, Bailey and Arno [6]. Letting
$r=\tfrac{1}{\sqrt{2}}$, $a=3$ and $b=8$, our theorem states that

$\frac{\pi}{4}=\arg\left(1+\tfrac{1}{\sqrt{2}}e^{\frac{3}{8}2\pi
i}\right)=\operatorname{BBP}(1,16,8,A)=\sum_{k=0}^{\infty}\frac{1}{16^{k}}\sum_{j=1}^{8}\frac{a_{j}}{8k+j},$

where
$A=(a_{1},\dots,a_{8})=(\frac{1}{2},\frac{1}{2},\frac{1}{4},0,\frac{-1}{8},\frac{-1}{8},\frac{-1}{16},0)$.

We can interpret the main result geometrically in the following way. Consider
the circle in $\mathbb{C}$ with centre $1$ and radius $r$, along with the
triangle subtended by the three points $A=0$, $B=1$ and $C=1+re^{i\theta}$.
Clearly, $C$ lies on the circle considered. Moreover, the angle at $B$ is
$\pi-2\theta\pi$, a rational multiple of $\pi$. If also the angle at $A$ is a
rational multiple of $\pi$, the theorem yields a BBP-type formula for $\pi$.
In the example given, the triangle is an isosceles triangle with angles at $A$
and $B$ equal to $\frac{\pi}{4}$, see figure 1.

$\tfrac{\pi}{4}$$\tfrac{\pi}{4}$$A=0$$C=1+\tfrac{1}{\sqrt{2}}e^{\frac{3}{8}2\pi
i}$$B=1$

Figure 1. Proof that
$\frac{\pi}{4}=\operatorname{BBP}(1,16,8,(\frac{1}{2},\frac{1}{2},\frac{1}{4},0,\frac{-1}{8},\frac{-1}{8},\frac{-1}{16},0))$

Of course, such triangles are rather uncommon. However, as we will see, we can
sometimes combine triangles, which are not necessarily as nice, and obtain new
BBP-type formulas. The observation which makes this work is the following
lemma.

######  Lemma 2.

For fixed $d,n\in\mathbb{N}$ and $b>1$, the map
$\operatorname{BBP}(d,b,n,\cdot):\mathbb{R}^{n}\rightarrow\mathbb{R}$ is
linear.

The proof is simply a matter of manipulating convergent series and hence
ommitted.

We would like to consider all BBP-type formulas arising from the main theorem
and look for combinations of them, which give formulas of interest. Recall
that the formulas arising in this way are for the angle of the $ABC$-triangle
at $A$, so unless this number is of interest, the formulas obtained are not
interesting in themselves, although it may happen that they combine with other
formulas to reveal a formula of interest. In order to study these
combinations, we define a function which given an even natural number $b$ we
define the Circle-To-BBP map
$\operatorname{CTB}_{b}:(0,1)\times\mathbb{N}\rightarrow\mathbb{R}^{b}$ by
letting

$\operatorname{CTB}_{b}(r,a)=A=(a_{1},\dots,a_{b}),\quad\text{with
}a_{j}=r^{j}(-1)^{j+1}\sin\left(\tfrac{ja}{b}2\pi\right).$

The map reads off the vector associated to the BBP-type formula arising from
Theorem 1 with the parameters $r$, $b$ and $\theta=\frac{a}{b}2\pi$. In other
words, the range of the map is the collection of permissible vectors from
triangles arising by letting the point $C$ vary over all corners of a regular
$b$-gon inscribed in a circle with centre $1$ and radius $r$.

By periodicity of sine, $\operatorname{CTB}_{b}(r,\cdot)$ is clearly
$b$-periodic, but in fact one only needs to consider
$a\in\\{1,\dots,\frac{b}{2}-1\\}$. Indeed, using that sine is odd and
$2\pi$-periodic, we immediately verify the relations in the following lemma.

######  Lemma 3.

For $b$ an even natural number and $r\in(0,1)$,

  1. (1)

$\operatorname{CTB}_{b}(r,0)=\operatorname{CTB}_{b}(r,\frac{b}{2})=0$.

  2. (2)

$\operatorname{CTB}_{b}(r,a)=-\operatorname{CTB}_{b}(r,b-a)$.

##  3\. An algorithm

At this point, we are ready to turn the above results into an algorithm. We
first explain how our algorithm works when searching for formulas for $0$. We
will subsequently explain how to modify it to look for formulas for other
constants. In this case, algorithm works as follows.

  1. (1)

Fix an algebraic base number $r\in(0,1)$ (strictly speaking, the base will be
$r^{-1}$) and an even natural number $b$.

  2. (2)

Compute all vectors $\operatorname{CTB}_{b}(r,a)$ for
$a\in\\{1,\dots,\frac{b}{2}-1\\}$.

  3. (3)

Using a standard procedure [8] for finding approximate integer relations
between these vectors. We have used the LLL implementation in sage in our
algorithm. We describe how to do this below.

  4. (4)

Verify that any approximate identities are in fact proper identities. In our
case, we have used exact algebraic arithmetic in sage for this purpose on
noting that all numbers in the vectors $\operatorname{CTB}_{b}(r,a)$ have
coordinates which are arguments of algebraic numbers.

  5. (5)

Output all relations found in this way.

  6. (6)

Check for linear dependence among the obtained formulas to ensure that they
are properly distinct. This is done by checking whether the Gram matrix of any
selection of vectors is zero.

The procedure in [8] describes how to find approximate integer relations
between a set of real numbers. What we really need is an approximate relation
between vectors, but in our case, this is the same thing if we work instead
with the BBP-type formulas. Indeed, we do the following. First, we let
$x_{i}=\operatorname{BBP}(1,r^{-b},b,\operatorname{CTB}_{b}(r,i))$ for
$i=1,\dots,\frac{b}{2}-1$, and take suitably good rational approximations
$\hat{x}_{i}$ to these numbers. We now pick a large integer $N$ and let
$m_{i}=[N\hat{x}_{i}]$, the nearest integer to $N\hat{x}_{i}$. Consider now
the lattice $\Lambda\subseteq\mathbb{R}^{\frac{b}{2}}$, spanned by the vectors
$v_{i}=(0,\dots,0,1,0,\dots,m_{i})$ with $1$ in the $i$’th coordinate and
$m_{i}$ in the $n+1$’st coordinate. Applying LLL to this lattice yields a non-
zero vector of small norm. But any vector of the lattice $\Lambda$ is of the
form

$\sum_{i=1}^{b/2-1}a_{i}v_{i}=\left(a_{1},\dots,a_{b/2-1},\sum_{i=1}^{b/2-1}a_{i}m_{i}\right).$

This vector is short only if the final coordinate is small, that is, if we
have a near integer relation between the $m_{i}$.

We would hope that the relation is in fact an integer relation between the
$x_{i}$, i.e. that $\sum_{i=1}^{b/2-1}a_{i}x_{i}=0$. This is not an
unreasonable hope. Indeed, if we let $\epsilon>0$ be the worst among the
qualities of our approximations $\hat{x}_{i}$, then

$\displaystyle N\left|\sum_{i=1}^{b/2-1}a_{i}x_{i}\right|$
$\displaystyle=\left|\sum_{i=1}^{b/2-1}a_{i}N(x_{i}-\hat{x}_{i}+\hat{x}_{i})\right|$
$\displaystyle\leq
N\epsilon\left|\sum_{i=1}^{b/2-1}a_{i}\right|+\left|\sum_{i=1}^{b/2-1}a_{i}N\hat{x}_{i}\right|$
$\displaystyle\leq
N\epsilon\left|\sum_{i=1}^{b/2-1}a_{i}\right|+\left|\sum_{i=1}^{b/2-1}a_{i}m_{i}\right|+\frac{\left|\sum_{i=1}^{b/2-1}a_{i}\right|}{2}.$

Letting

$s=\left\|\sum_{i=1}^{b/2-1}a_{i}v_{i}\right\|^{2}$

the square of the Euclidean norm of the vector found by LLL, this immediately
implies that

$\left|\sum_{i=1}^{b/2-1}a_{i}x_{i}\right|\leq\frac{s}{N}+|s|\epsilon.$

Note that by [8], the number $s$ is within a constant factor of the square of
the smallest non-zero Euclidean norm occurring for a vector in $\Lambda$. In
other words, for $\epsilon$ small and $N$ large, we will expect the the left
hand side to be really small and hence equal to zero to within some tolerance.
Whether this approximate equality is in fact an equality is subsequently
checked in the next step of the algorithm.

All that remains in order to include formulas for other numbers than zero is
to include them among the $x_{i}$. If one is interested in obtaining a formula
for $\pi$, say, one would just include $\pi$ in the list of $x_{i}$ as the
$b/2$’th element, so that the lattice $\Lambda$ will now lie in
$\mathbb{R}^{b/2+1}$. This will allow the algorithm to look for formulas for
linear combinations including the new number as well, and if one was looking
for integer combinations of other constants, these could also be added. Of
course, it is entirely possible that the integer relations coming out of the
search algorithm do not include the constants with a non-zero coefficient, in
which case we have just recovered a zero formula. In fact, if one plugs in an
entirely arbitrary number, in all probability the number in question does not
admit a BBP-type formula at all, and thus will be ignored by the algorithm
(though its inclusion may increase run time). In our search, we have looked
for formulas for $\pi$ and multiples thereof. This is due to the argument
below.

We now finally describe how to reduce a candidate identity to a problem in
exact algebraic arithmetic, cf. Step 4 of the algorithm. As described here,
the procedure works for $\pi$.

A candidate identity obtained in Step 3 has the form

$\sum_{i=1}^{N}a_{i}\arg\left(q+re^{i\theta_{i}}\right)=b\pi,$

where the $a_{i}$ and $b$ are integers, not all zero; $N$ is the length of the
formula, we are searching for (in our case $b/2-1$ or $b/2$); $r$ is
algebraic; and the $\theta_{i}$ are rational multiples of $\pi$. Thus,
$c_{i}=q+re^{i\theta_{i}}\in\overline{\mathbb{Q}}$. Thus, rewriting and using
properties of the argument,

$b\pi=\sum_{i=1}^{N}a_{i}\arg(c_{i})=\sum_{i=1}^{N}\arg\left(c_{i}^{a_{i}}\right)=\arg\left(\prod_{i=1}^{N}c_{i}^{a_{i}}\right).$

Since Step 3 outputs an integer $b$, this implies that the number inside the
argument must be a real number, whence it must be equal to its complex
conjugate. The $a_{i}$ are likewise integers, so the equality to be verified
exactly is

$\prod_{i=1}^{N}c_{i}^{a_{i}}=\prod_{i=1}^{N}\overline{c_{i}}^{a_{i}},$

an identity in algebraic numbers, which can be verified in exact algebraic
arithmetic.

An implementation of the above algorithm in sage can be found on
https://projects.au.dk/ada/

The algorithm described above outputs any BBP-type zero formulas found in this
way. As we will see below, it does not recover all known formulas.
Nonetheless, we believe that the method has its merits. One such feature is
the easily understood geometrical interpretation of the algorithm. Another is
the fact that $r$ need not be an integer, but only an algebraic number. In our
experiments, we have attempted to look for formulas with a Pisot number base.

Another feature of the algorithm is that the format of the output in fact
constitutes a geometric proof of the formula in question. This is to be
understood in the sense that the output is a linear combination of BBP-type
formulas given in terms of the $\operatorname{CTB}$-function. In order to
remove any computer interaction from the proof of a given formula, one only
needs to verify that the associated triangles have a good configuration. This
is best illustrated with an example. The null formula

(2)
$\begin{split}0=2^{12}\big{(}&\operatorname{BBP}(1,2^{12},24,\operatorname{CTB}_{24}(1/\sqrt{2},5))\\\
&-\operatorname{BBP}(1,2^{12},24,\operatorname{CTB}_{24}(1/\sqrt{2},11))\big{)},\end{split}$

was originally found by Lafont in [7], though in the usual BBP-notation. Out
algorithm recovers this formula. The form of the formula above allows us to
turn the formula above into an easily verifiable geometric condition as
illustrated in Figure 2. The large triangle gives rise to the first term,
while the smaller triangle gives rise to the second term. The proof of the
formula is simply a matter of verifying the fact that the complex numbers $0$,
$1+2^{-1/2}e^{\frac{5}{12}\pi i}$ and $1+2^{-1/2}e^{\frac{11}{12}\pi i}$ are
collinear – essentially an exercise for undergraduates.

$A=0$$C=1+\tfrac{1}{\sqrt{2}}e^{\frac{5}{24}2\pi
i}$$C^{\prime}=1+\tfrac{1}{\sqrt{2}}e^{\frac{11}{24}2\pi i}$$B=1$

Figure 2. Geometric proof of Lafont’s null formula.

##  4\. Computational results

We consider first integer bases, which are powers of $2$ and $3$. In Table 1,
we give a formula count of the known formulas in various bases and lengths
with integer coefficients. In addition to the formulas with integer
coefficients, our algorithm also finds formulas for $0$ with non-integer
coefficients. These are omitted from the counts. All formulas obtained are of
degree $1$. We list the values of base, the associated $r$ and the length $b$
as they occur in Theorem 1.

base | $r$ | $b$ | # known formulas | # recovered formulas
---|---|---|---|---
$2^{4}$ | $1/\sqrt{2}$ | $8$ | $1$ | $0$
$2^{6}$ | $1/2$ | $6$ | $1$ | $0$
$2^{12}$ | $1/\sqrt{2}$ | $24$ | $5$ | $1$
$2^{20}$ | $1/\sqrt{2}$ | $40$ | $3$ | $1$
$3^{6}$ | $1/\sqrt{3}$ | $12$ | $2$ | $1$
Table 1. Null formulas in integer bases

In our condensed notation using the $\operatorname{CTB}$-notation which allows
for easy geometric proof of the formulas, we have obtained the following known
formulas aside from Lafont’s formula (2) discussed above.

$\begin{split}0=2^{20}\big{(}&\operatorname{BBP}(1,2^{20},40,\operatorname{CTB}_{40}(1/\sqrt{2},3))\\\
&+\operatorname{BBP}(1,2^{20},40,\operatorname{CTB}_{40}(1/\sqrt{2},5))\\\
&+\operatorname{BBP}(1,2^{20},40,\operatorname{CTB}_{40}(1/\sqrt{2},11))\\\
&-\operatorname{BBP}(1,2^{20},40,\operatorname{CTB}_{40}(1/\sqrt{2},13))\\\
&-\operatorname{BBP}(1,2^{20},40,\operatorname{CTB}_{40}(1/\sqrt{2},15))\\\
&+\operatorname{BBP}(1,2^{20},40,\operatorname{CTB}_{40}(1/\sqrt{2},19))\\\
&-\operatorname{BBP}(1,2^{20},40,\operatorname{CTB}_{40}(1/\sqrt{2},20))\big{)},\end{split}$

which one finds in [3]; and

(3)
$\begin{split}0=386\sqrt{3}\big{(}&\operatorname{BBP}(1,3^{6},12,\operatorname{CTB}_{12}(1/\sqrt{3},3))\\\
&-\operatorname{BBP}(1,3^{6},12,\operatorname{CTB}_{12}(1/\sqrt{3},5))\big{)},\end{split}$

found by Adegoke in [1]. The leading coefficients, which admittedly look a
little silly, are there to ensure that if one rewrites the obtained formula as
a single BBP-type formula, it becomes a pure formula without leading
coefficient and with integer coefficients. For instance, on expanding and
regrouping terms, (3) becomes

$0=\operatorname{BBP}(1,729,12,(243,-243,-324,-81,27,0,-9,9,12,3,-1,0)).$

Using our algorithm, we have found a BBP-type formula for $\sqrt{3}\pi$, which
we have not been able to find in litterature. The formula reads

$\sqrt{3}\pi=\frac{9}{2^{12}}\operatorname{BBP}(1,2^{12},24,A),$

with

$\begin{split}A=(&2^{11},0,0,2^{10},2^{9},0,2^{8},2^{8},0,0,2^{6},0,\\\
&-2^{5},0,0,-2^{4},-2^{3},0,-2^{2},-2^{2},0,0,-1,0).\end{split}$

The computer generated proof found by the algorithm is

$\begin{split}\sqrt{3}\pi=3\sqrt{3}\big{(}&\operatorname{BBP}(1,2^{12},24,\operatorname{CTB}_{24}(1/\sqrt{2},5))\\\
&+\operatorname{BBP}(1,2^{12},24,\operatorname{CTB}_{24}(1/\sqrt{2},11))\big{)}.\end{split}$

The BBP-type formula is found by expanding the right hand side.

The real merit of our algorithm is that it is able to produce BBP-type
formulas to non-integer bases, while still producing formulas with integer
coefficients. We have not been able to find examples of this in literature,
though Adegoke [2] does have BBP-type formulas in base $\phi$, the Golden
ratio. However, the formulas obtained by Adegoke have coefficients in
$\mathbb{Z}[\phi]$, which to us seems a little unnatural. The bases studied
must be algebraic. We have made our algorithm run with several Pisot bases. Of
course, the output will always be expansions with rational coefficients as is
the case with integer bases, and as such they are not $\beta$-expansions in
the sense of Rényi [10] and Parry [9]. Nevertheless, we find the expansions of
interest. We are particularly intrigued by the rather large number of
expansions with the Golden Ratio as a base.

In Table 2, we present the number of null formulas found to various Pisot
bases of low degree and height. In each case, the minimal polynomial of the
Pisot number; the approximate value of the Pisot number, i.e. the real root
which is greater than one $1$, is given; the denominator $b$ is given; and the
number of null formulas found are all given. The value of $r$ is in each case
set equal to the reciprocal of the Pisot number in question. The values of $b$
chosen are $60$ and $24$ due to the smoothness of these numbers. The effect is
that our algorithm finds formulas coming from even divisors of the numbers.

Polynomial | Approximate value | $b$ | # formulas obtained
---|---|---|---
$x^{2}-x-1$ | $1.618033$ | $60$ | $12$
$x^{2}-2x-2$ | $2.732050$ | $60$ | $1$
$x^{2}-2x-1$ | $2.414213$ | $60$ | $0$
$x^{2}-3x+1$ | $2.618033$ | $60$ | $4$
$x^{2}-3x-1$ | $3.302775$ | $60$ | $0$
$x^{2}-4x+2$ | $3.414213$ | $60$ | $0$
$x^{3}-x-1$ | $1.3247179$ | $60$ | $3$
$x^{4}-x^{3}-1$ | $1.380277569$ | $24$ | $1$
Table 2. Null formulas in Pisot bases

The first two lines of Table 2 are the Golden Ratio and $\sqrt{3}+1$
respectively. We denoting the Golden ratio by $\phi$, and putting
$\psi=\sqrt{3}+1$, we present a small selection of the formulas obtained by
our algorithm. For $\phi$, we have for instance the following formulas. Since
searching for BBP-type formulas in non-integer bases appear to be largely
unexplored territory, we have not found any of these in literature.

$\begin{split}0=-\big{(}&\operatorname{BBP}(1,\phi^{60},60,\operatorname{CTB}_{60}(\phi^{-1},18))\\\
&+\operatorname{BBP}(1,\phi^{60},60,\operatorname{CTB}_{60}(\phi^{-1},24)),\end{split}$
$\begin{split}0=-\big{(}&\operatorname{BBP}(1,\phi^{60},60,\operatorname{CTB}_{60}(\phi^{-1},14))\\\
&+\operatorname{BBP}(1,\phi^{60},60,\operatorname{CTB}_{60}(\phi^{-1},26)),\end{split}$
$\begin{split}0=-\big{(}&\operatorname{BBP}(1,\phi^{60},60,\operatorname{CTB}_{60}(\phi^{-1},8))\\\
&+\operatorname{BBP}(1,\phi^{60},60,\operatorname{CTB}_{60}(\phi^{-1},28)),\end{split}$

and

$\begin{split}0=-\big{(}&\operatorname{BBP}(1,\phi^{60},60,\operatorname{CTB}_{60}(\phi^{-1},8))\\\
&-\operatorname{BBP}(1,\phi^{60},60,\operatorname{CTB}_{60}(\phi^{-1},9))\\\
&+\operatorname{BBP}(1,\phi^{60},60,\operatorname{CTB}_{60}(\phi^{-1},21)).\end{split}$

In addition, our algorithm found a formula for for $\pi$ to base $\phi$,

$\pi=10\operatorname{BBP}(1,\phi^{15},15,\operatorname{CTB}(\phi^{-1},2)).$

Note that the $b$ in the latter formula is reduced to $15$. This reduction of
the base is possible by the above remarks on the smoothness of $60$.

For $\psi=\sqrt{3}+1$, the following null formula occurs in base $12$, even
though we performed our search in base $60$,

$\begin{split}0=\big{(}&\operatorname{BBP}(1,\psi^{12},12,\operatorname{CTB}_{12}(\psi^{-1},2))\\\
&-\operatorname{BBP}(1,\psi^{12},12,\operatorname{CTB}_{12}(\psi^{-1},5)).\end{split}$

While this is the only null formula, which we have found, a formula for $\pi$
did show up,

$\pi=12\operatorname{BBP}(1,\psi^{12},12,\operatorname{CTB}_{12}(\psi^{-1},2)).$

Of course, by the two preceding formulas, it immediately follows that

$\pi=12\operatorname{BBP}(1,\psi^{12},12,\operatorname{CTB}_{12}(\psi^{-1},5)).$

As a final remark on the non-integer formulas obtained, for the number
$\phi+1$ with minimal polynomial $x^{2}-3x+1$, our usual approximation
parameter, measured by $N$ in the application of LLL, of $N=10^{10}$ only
resulted in two null formulas. Increasing the parameter to $N=10^{15}$
resulted in an additional two formulas. This suggests that we may be missing
formulas, which should nevertheless show up as arising from triangles in the
way described in this paper, due to the restrictions in sage.

##  5\. Concluding remarks

We end this paper with some concluding remarks and open problems.

First, the abundance of zero formulas found for the Golden Ratio is somewhat
of a mystery! We really do not know why this particular irrational base yields
so many zero formulas. Judging from the other Pisot numbers studied, it
appears not to have anything to do with the Pisot’ness of the number. Of
course, it is entirely possible that we have just not found all zero formulas
for other numbers. The issue definitely requires further study.

Second, the original intention behind BBP-type formulas was to find deep
digits in the expansions of numbers such as $\pi$, and indeed the original BBP
formula allows one to calculate the $n$’th binary digit of $\pi$ without
knowing all the preceding ones. For irrational bases, such as the Pisot bases
studied in this paper, one also has a digital expansion, namely the
$\beta$-expansions of Rényi [10] and Parry [9]. It is not clear to us how to
extract digital information from the expansions produced by our setup. Naively
trying to extend the spigot algorithms to the irrational setup fails
miserably. Obtaining digital information from our formulas remain an open
problem.

Third and finally, we are only able to deal with BBP-type formulas of degree
one in the present setup. Obtaining a similarly intuitive and geometrical
approach to formulas of higher degree also remains an open problem. It would
appear that the simple interpretation with triangles with vertices at $0$, $1$
and a point on a regular polygon on a particular circle is no longer
appropriate for this type of problem. We leave it for the interested reader to
ponder how to extend the interpretation of the formulas and their
algorithmical implementation to the higher degree case.

## References

  * [1] Kunle Adegoke.  New binary and ternary digit extraction (BBP-type) formulas for trilogarithm constants.  New York J. Math., 16:361–367, 2010.
  * [2] Kunle Adegoke.  The golden ratio, Fibonacci numbers and BBP-type formulas.  Fibonacci Quart., 52(2):129–138, 2014.
  * [3] D. Bailey.  A compendium of BBP-type formulas for mathematical constants.  `https://www.davidhbailey.com/dhbpapers/bbp-formulas.pdf\verb`, 2017\.
  * [4] David Bailey, Peter Borwein, and Simon Plouffe.  On the rapid computation of various polylogarithmic constants.  Math. Comp., 66(218):903–913, 1997.
  * [5] Daniel Barsky, Vicente Muñoz, and Ricardo Pérez-Marco.  On the genesis of BBP formulas.  Acta Arith., 198(4):401–426, 2021.
  * [6] Helaman R. P. Ferguson, David H. Bailey, and Steve Arno.  Analysis of PSLQ, an integer relation finding algorithm.  Math. Comp., 68(225):351–369, 1999.
  * [7] J. O. Lafont.  Obtention of a zero relation of the form $0={P}(1,2^{12},24,{A})$.  `http://oeis.org/wiki/User:Jaume_Oliver_Lafont/BBP`.
  * [8] A. K. Lenstra, H. W. Lenstra, Jr., and L. Lovász.  Factoring polynomials with rational coefficients.  Math. Ann., 261(4):515–534, 1982.
  * [9] W. Parry.  On the $\beta$-expansions of real numbers.  Acta Math. Acad. Sci. Hungar., 11:401–416, 1960.
  * [10] A. Rényi.  Representations for real numbers and their ergodic properties.  Acta Math. Acad. Sci. Hungar., 8:477–493, 1957.


# Probing the Viscoelastic Properties of Aqueous Protein Solutions using
Molecular Dynamics Simulations

D. F. Hanlon<EMAIL_ADDRESS>Ivan Saika-Voivod  G. Todd Andrews  Department of
Physics and Physical Oceanography, Memorial University of Newfoundland, St.
John’s, NL, Canada, A1B 3X7  M. Shajahan G. Razul  Department of Chemistry,
St. Francis Xavier University, Antigonish, NS, B2G 2W5, Canada

###### Abstract

We performed molecular dynamics simulations to investigate the viscoelastic
properties of aqueous protein solutions containing an antifreeze protein, a
toxin protein, and bovine serum albumin. These simulations covered a
temperature range from 280 K to 340 K. Our findings demonstrate that lower
temperatures are associated with higher viscosity as well as a lower bulk
modulus and speed of sound for all the systems studied. Furthermore, we
observe an increase in the bulk modulus and speed of sound as the temperature
increases up to a weak maximum while the viscosity decreases. Moreover, we
analyzed the influence of protein concentration on the viscoelastic properties
of the antifreeze protein solution. We observed a consistent increase in the
bulk modulus, speed of sound, and viscosity as the protein concentration
increased. Remarkably, our molecular dynamics simulations results closely
resemble the trends observed in Brillouin scattering experiments on aqueous
protein solutions. The similarity thus validates the use of simulations in
studying the viscoelastic properties of protein water solutions. Ultimately,
this work provides motivation for the integration of computer simulations with
experimental data and holds potential for advancing our understanding of both
simple and complex systems.

###### pacs:

Valid PACS appear here

††preprint: APS/123-QED

##  I Introduction

Water is vital for life on Earth yet there are still many properties of water
and its interaction with other molecules that are not fully understood. In
particular, there has been considerable recent interest in the viscoelastic
properties of aqueous biomacromolecular solutions, especially water-protein
solutions, due to their importance in understanding fundamental biophysical
and biochemical processes.

Various experimental techniques have been used to study the viscoelastic
properties of aqueous protein solutions Hanlon  _et al._ (2023a); Denny
(1984); Sharma  _et al._ (2011); Dmitriev  _et al._ (2019). A previous
Brillouin scattering study examined aqueous solutions of lysozyme and bovine
serum albumin (BSA) with varying concentrations Dmitriev  _et al._ (2019). The
study found that increasing protein concentration increases both sound
velocity and viscosity, determined by spectral peak position and peak full
width at half maximum (FWHM), respectively. In a rheology study, the bulk and
viscoelastic properties of BSA at low concentrations were investigated Sharma
_et al._ (2011). Both the Brillouin and rheology studies showed a strong
dependence of viscoelasticity on protein concentration. Similarly, another
study used Brillouin light scattering spectroscopy to examine the temperature
dependence of the viscoelastic properties of snail mucus, a natural system
consisting primarily of water and glycoprotein Hanlon  _et al._ (2023a). As
the temperature increased, the speed of sound and storage modulus increased,
while viscosity and hypersound attenuation decreased. Another notable feature
present in the mucus was a liquid-to-solid phase transition at around 270.5 K.
Further experiments on diluted and dehydrated snail mucus showed significant
variations in speed of sound, viscosity, and sound absorption as a function of
protein concentration Hanlon  _et al._ (2023b). Rheology studies were also
conducted on dehydrating snail mucus and revealed that the shear modulus
increased with increasing degree of dehydration Denny (1984). Overall,
previous experimental studies have consistently demonstrated that the
viscoelasticity of aqueous protein solutions is strongly influenced by protein
concentration Hanlon  _et al._ (2023a); Denny (1984); Sharma  _et al._ (2011);
Dmitriev  _et al._ (2019).

Somewhat surprisingly, molecular dynamics (MD) simulations have not been used
much for investigating bulk viscoelastic properties such as bulk modulus,
speed of sound, and viscosity of aqueous protein solutions nor fluids in
general. Only a few MD studies have explored the dynamical properties of
aqueous protein solutions containing collagen molecules and different co-
solvents Gautieri  _et al._ (2012); Wei  _et al._ (2010). Specifically, MD
simulations have been employed to determine the viscosity and relaxation times
of collagen molecules, aiming to explore their elastic behavior and
investigate the possible role of cross-linked collagen molecules in
viscoelastic properties, as previously hypothesized Gautieri  _et al._ (2012);
Silver  _et al._ (2001a, b). A likely reason for the limited studies on the
viscoelastic properties of aqueous protein solutions is possibly due to the
overall accuracy and reliability of the simulations which depend on a number
of factors, including the choice of force field, system size, and simulation
time.

In this paper we demonstrate the capabilities of MD simulations in exploring
fundamental properties of fluids by obtaining bulk viscoelastic properties of
protein-water systems over the temperature range 280 K $\leq T\leq$ 340 K. We
found that the magnitude and temperature dependence of the bulk modulus, sound
speed, and shear viscosity of simulated aqueous protein solutions containing
an antifreeze protein, a toxin protein, and BSA are similar to those observed
in previous experimental studies on aqueous protein solutions. The combination
of MD simulations with experimental results therefore provides a more
comprehensive understanding of the physics of molecular interaction in water-
macromolecule solutions that could lead to advancements in fields such as drug
discovery, materials science, biochemistry and biophysics.

##  II Computational Details

Table 1 gives an overview of the simulated systems studied in the present
work. These included aqueous solutions of an antifreeze glycoprotein found in
eelpout, a cross-linking toxin protein found in Vibrio cholerae, and the
widely studied bovine serum albumin Jia  _et al._ (1996); Durand  _et al._
(2012); Bujacz (2012). The initial protein configurations were obtained from
the Protein Data Bank (PDB). The naming convention employed to designate each
system consists of the PDB ID of the protein followed by the weight percentage
of that protein in the system. For instance, antifreeze protein 1MSI at 4 wt%
is referred to as 1MSI-4. The antifreeze glycoprotein 1MSI was used due to its
simplistic structure and the fact that antifreeze glycoprotein solutions have
recently been the subject of experimental work Hanlon  _et al._ (2023b, a).
The toxin protein with PDB ID 4GQK was used primarily for its ability to
cross-link. Lastly, BSA (PDB ID 4F5S) was used since it is an extensively
studied molecule. Moreover, in addition to the specific reasons for using the
chosen proteins, their collective inclusion enables a comprehensive study of
proteins characterized by progressively higher molecular weights. This allows
for the potential to provide insights into the influence of protein size on
the viscoelastic properties of such systems.

Table 1: Overview of protein-water systems used in the present molecular dynamics simulations. The temperature range studied was 280 K - 340 K. System $\rightarrow$ | 1MSI/Water111System naming convention: Protein Data Bank ID followed by protein concentration in wt% (e.g., 1MSI-4 is protein 1MSI at 4 wt%). | 4GQK/Water | 4F5S/Water
---|---|---|---
Protein Data Bank ID | 1MSI | 4GQK | 4F5S
Protein |  Antifreeze Glycoprotein QAE(HPLC 12) |  VgrG1-ACD with ADP |  Bovine Serum Albumin
Weight of Protein (kDa) | 7.4 | 44.3 | 133.3
# of Protein Molecules in System | 1, 3, 5 | 2 | 1
# of Water Molecules in System | 9970, 8110, 7190 | 52,030 | 258,840
Protein Concentration (wt %) | 4, 12, 22 | 4 | 3
System Volume (nm3) | 6.5$\times$6.5$\times$6.5 | 12$\times$12$\times$12 | 20$\times$20$\times$20
System Density222Density of each system at 300 K. (g/cm3) | 1.009, 1.025, 1.060 | 1.024 | 1.010
Rationale for Use |  Simple structure & an antifreeze glycoprotein |  Ability to cross-link in solution |  Availability & a widely studied protein
Miscellaneous Notes |  Antifreeze glycoprotein found in eelpout Jia  _et al._ (1996) |  Toxin protein found in Vibrio cholerae Durand  _et al._ (2012) |  Structure determined from Ref. Bujacz (2012)

The MD simulations were carried out in the isothermal-isobaric (NPT) ensemble
using GROMACS v2022.3 Hess  _et al._ (2008); Lindahl  _et al._ (2001); Van Der
Spoel  _et al._ (2005). We used the TIP4P/2005 water model (volume $4\times
4\times 4$ nm3, density 0.992 g/cm3) as a control system due to its ability to
accurately represent the properties of water over the temperature range
studied Abascal and Vega (2005); González  _et al._ (2016); Yu  _et al._
(2023). Proteins in the system interact using the All-Atom Optimized Potential
for Liquid Simulations (OPLS-AA) Jorgensen  _et al._ (1996) and was chosen due
to its good description of liquid organic systems and ability to reproduce
thermodynamic properties such as density and heat of vaporization, structural
properties, as well as free energies of hydration Zangi (2018). The equations
of motion were integrated with a time step of 2 fs using the leap-frog
algorithm Cuendet and van Gunsteren (2007). Temperature was held constant
using a Nosé-Hoover thermostat with a time constant of 0.1 ps. The pressure
was held constant using the Parrinello-Rahman barostat with a time constant of
2 ps. The viscoelastic properties of each system (see Table 1) were
investigated over the temperature range 280 K to 340 K. For each system, the
size of simulation box and number of solute molecules vary due to the
differing sizes of the proteins used and the desire to replicate approximate
protein concentrations used in experimental work.

##  III Results & Discussion

###  III.1 Bulk Modulus

Figure 2 shows the bulk modulus ($K_{s}$) for solutions studied in this work
(see Table 1) obtained by analyzing the fluctuations in volume and pressure in
the simulations. For most systems, the bulk modulus increases with increasing
temperature until it reaches a maximum value near 330 K. However, exceptions
to this trend are observed in 1MSI-12 and 1MSI-22, where the bulk modulus
exhibits a maximum at 320 K and 300 K, respectively. The difference in the
maximum for these two solutions is caused by the increased protein
concentration in each system. Additionally, in the case of 1MSI solutions,
higher protein concentrations lead to a larger bulk modulus. The presence of a
maximum in the bulk modulus also suggests that the molecular structure in each
system is most rigid at this temperature.

Figure 1: Representative simulation box of the 1MSI-4 system containing a
single 1MSI protein (at center of box) and 2417 water molecules (grey flecks).
Figure was produced using VMD v.1.9 Humphrey  _et al._ (1996).  Figure 2:
Temperature dependence of bulk modulus for water and aqueous protein solutions
with different protein species and concentration. Figure 3: Temperature at
which maximum bulk modulus occurs (T${}^{*}_{K_{s}}$) as a function of
concentration for 1MSI solutions. The line of best fit, represented by the
dashed line with the equation T${}^{*}_{K_{s}}$ = 337[Xp] - 1.68, where [Xp]
denotes the protein concentration.

Figure 3 illustrates the temperature at which the bulk modulus is at a maximum
for each 1MSI solution. It is evident that as protein concentration increases,
the temperature at which the maximum bulk modulus occurs decreases. The reason
behind this phenomenon is due to the fact that the behavior of water and its
bulk modulus (compressibility) is influenced by hydrogen bonds between water
molecules and proteins and the presence of proteins has been known to disrupt
these bonds Ball (2008).

###  III.2 Sound Velocity

Figure 4: Temperature dependence of sound velocity of water and aqueous
protein solutions with different protein species and concentration.

The temperature dependence of the sound velocity for each of the protein-water
systems was determined directly from the bulk modulus and the solution mass
density $\rho$ using the relation

$v=\sqrt{\frac{K_{s}}{\rho}}.$ (1)

As shown in Fig. 4, the behaviour is qualitatively similar for all systems.
There is a steady increase in sound velocity with increasing temperature for
temperatures at the lower end of the investigated range while at higher
temperatures the velocity decreases with increasing temperature. The maximum
in sound velocity occurs at $\sim 330$ K for systems with low protein
concentrations (3-4%) and at slightly lower temperatures for those with higher
concentrations (12% and 22%).

In the case of the 1MSI systems, the sound velocity displays a strong
dependence on concentration. Specifically, an increase in protein
concentration generally leads to a higher sound velocity, except for the
values at 340 K, where 1MSI-4 exhibits a larger sound velocity than 1MSI-12.
The reason for this unexpected behavior at 340 K is currently unknown.
Additionally, as the concentration increases, the maximum sound velocity for
the 1MSI systems occurs at lower temperatures. This can be seen in Fig. 3
which shows the temperature at which the maximum bulk modulus occurs, directly
related to the speed of sound via Eq. 1. This indicates a complex interplay
between concentration, temperature, and sound velocity within the 1MSI system.
Notably, except at the highest temperatures, the sound velocities for the 1MSI
systems are larger than those for water. Conversely, the 4F5S-3 and 4GQK-4
systems exhibit sound velocity values that closely resemble those of water,
although these velocities are generally lower than those of water, despite
having protein concentrations similar to 1MSI-4. This observation is
interesting because it demonstrates the impact of protein species on sound
velocity in aqueous protein solutions, ranging from a decrease relative to
v${}_{H_{2}O}$ to an increase relative to v${}_{H_{2}O}$.

The sound velocities for 1MSI-4 and 1MSI-12 systems are similar for all
temperatures but quite different from those of the 1MSI-22 system.
Furthermore, the sound velocity for the 1MSI-22 system is also much greater
than for all other systems, the difference being most obvious at lower
temperatures ($T\leq$ 300 K). This overall larger sound velocity for the
1MSI-22 system can be understood by noting that sound velocity is directly
related to bulk modulus. Therefore, with increased protein concentration, the
liquid becomes less compressible than at lower concentrations resulting in an
increased bulk modulus and consequently, sound velocity.

Another notable feature in the bulk modulus and sound velocity data is the
near convergence of values for all systems near 340 K (see Fig. 4). At higher
temperature this type of feature has been associated with protein denaturation
and is consistent with nuclear magnetic resonance studies and MD simulations
on lysozyme-water solutions that provide evidence for denaturation at $\sim
340$ K Mallamace  _et al._ (2011); Zhang  _et al._ (2009). However, its worth
noting that our investigation into the temperature dependent changes in the
radius of gyration did not reveal any signs of protein denaturation occurring.

###  III.3 Shear Viscosity

Figure 5: Shear viscosity of water and aqueous protein solutions versus
temperature for proteins 4F5S, 4GQK, and 1MSI at concentrations of 4%, 12%,
and 22%. Also shown for the purposes of comparison are apparent viscosity
values for water and snail mucus, a natural system consisting primarily of
water and glycoproteins, obtained from Brillouin scattering experiments Hanlon
_et al._ (2023a).

Figure 5 shows the shear viscosity as a function of temperature, calculated
using the Einstein relation Hess  _et al._ (2008)

$\eta=\lim_{t\to\infty}\frac{V}{k_{B}T}\frac{d}{dt}\Big{\langle}\Big{(}\int_{t_{0}}^{t_{0}+t}P_{xz}(t^{\prime})dt^{\prime}\Big{)}^{2}\Big{\rangle}_{t_{0}},$
(2)

where $V$ and $T$ are the volume and temperature of the system, respectively.
$P_{xz}(t^{\prime})$ is the off-diagonal component of the pressure tensor at
time $t^{\prime}$, $k_{B}$ is the Boltzmann constant, $t_{0}$ is the starting
time, and $t$ represents the integration time. Our results for the TIP4P/2005
water model show a shear viscosity value of 0.85 $\pm$ 0.09 cP at 300 K, which
is in good agreement with the previously obtained value of 0.855 cP at 298 K
from Ref. González and Abascal (2010).

Within uncertainty the viscosity of all systems decrease with increasing
temperature over the entire temperature range studied. This trend aligns with
experimental findings on similar protein-water systems Hanlon  _et al._
(2023a); Lupi  _et al._ (2011); Comez  _et al._ (2012); Pochylski  _et al._
(2005) (see apparent viscosity results for snail mucus, a natural water-
glycoprotein system in Fig. 5). For the 1MSI solutions, there is an non-linear
increase in shear viscosity as the temperature decreases. It is also clear
that the shear viscosity of the 1MSI solutions is strongly influenced by
protein concentration. In contrast, the magnitude and behaviour of the shear
viscosity of the simulated 4GQK-4 and 4F5S-3 solutions are similar to those of
TIP4P/2005.

Figure 6: Plot of natural logarithm of shear viscosity of simulated systems as a function of temperature. Solid lines indicate best fits for the entire temperature region. Also shown for the purposes of comparison is the natural logarithm of apparent viscosity for water and snail mucus, a natural system consisting primarily of water and glycoproteins, obtained from Brillouin scattering experiments Hanlon  _et al._ (2023a). Table 2: Best-fit parameters for fit of function $\displaystyle\ln(\eta)=\ln\eta_{0}+E_{a}/k_{B}T$ to computationally determined shear viscosity. System | $\ln\eta_{0}$ | $E_{a}$ | $R^{2}$
---|---|---|---
| (cP) | kJ/mol |
1MSI-4 | -5.9 $\pm$ 0.8 | 16.7 $\pm$ 0.7 | 0.987
1MSI-12 | -4.3 $\pm$ 0.9 | 14.5 $\pm$ 0.9 | 0.990
1MSI-22 | -11.5 $\pm$ 0.9 | 31.3 $\pm$ 0.8 | 0.977
4GQK-4 | -6.3 $\pm$ 0.6 | 18.2 $\pm$ 0.7 | 0.985
4F5S-3 | -3.6 $\pm$ 0.9 | 18.5 $\pm$ 0.9 | 0.977
TIP4P/2005 | -7.4 $\pm$ 0.6 | 18.5 $\pm$ 0.5 | 0.983
Water111Water previously studied by Ref. Hanlon  _et al._ (2023a) | -5.6 $\pm$ 0.3 | 13.8 $\pm$ 0.3 | 0.991
Mucus222Mucus previously studied by Ref. Hanlon  _et al._ (2023a) | -4.3 $\pm$ 0.3 | 13.6 $\pm$ 0.4 | 0.988
Water333Water previously studied by Ref. Lupi  _et al._ (2011) | — | 13.6 $\pm$ 0.6 | —

Figure 6 displays the natural logarithm of shear viscosity as a function of
inverse temperature for all systems in this study. All simulated systems show
a linear dependence on $1/T$ over the entire range studied. We therefore fit
the Arrhenius relationship

$\eta=\eta_{0}e^{E_{a}/k_{B}T}$ (3)

to the simulation data to obtain the activation energy $E_{a}$ and entropic
pre-factor $\eta_{0}$ (see Table 2). Notably, with the exception of 1MSI-22
the protein-water systems investigated in the current study have similar
values of activation energy, ranging from 16.7 to 18.5 kJ/mol. These findings
are consistent with previous experiments on aqueous protein and polymer
systems Lupi  _et al._ (2011); Comez  _et al._ (2012); Hanlon  _et al._
(2023a); Pochylski  _et al._ (2005), where the $E_{a}$ values were
approximately 14-24% lower than those observed in our simulations. To address
this difference we normalized the activation energy relative to ultra-pure
water for experimental data and the TIP4P/2005 water model for simulations.
This normalization allowed us to mitigate the influence of differences between
the ultra-pure water used in experiments and the water model used in
simulations. Figure 7 presents the ratios $E^{P}_{a}/E^{W}_{a}$ for the
protein-water solutions studied in this work and those obtained from Brillouin
spectroscopy for various dilutions of snail mucus Hanlon  _et al._ (2023b),
$E^{M}_{a}/E^{W}_{a}$. Remarkably, except for 1MSI-22, the ratios for the
protein-water systems analyzed in our research were close to unity. This
finding aligns with the behavior observed in natural snail mucus, where the
activation energy exhibited negligible variation across different
concentrations. It is also consistent with the results of a Brillouin
scattering study of water-tert-butyl alcohol in which activation energy was
found to be nearly independent of polymer concentration Lupi  _et al._ (2011).

Figure 7: Top Panel: Ratio of activation energy for the simulated aqueous
protein solution ($E_{a}^{P}$) to activation energy obtained for the
TIP4P/2005 water model ($E_{a}^{W}$). Bottom Panel: Ratio of activation energy
for hydrated snail mucus ($E_{a}^{M}$) to activation energy for ultra-pure
water ($E_{a}^{W}$), as reported in Ref. Hanlon  _et al._ (2023b). The labels
on the horizontal axis in the lower panel are approximate glycoprotein weight
percentages. Error bars are approximately the size of the symbols.

As evident from the results displayed in both Fig. 5 and 6, the simulated
systems show good agreement with both water and natural snail mucus as
previously obtained experimentally Hanlon  _et al._ (2023a). For example,
considering the 1MSI solutions, which contain antifreeze proteins of varying
concentrations, we see very similar behaviour in the viscosity of simulated
protein systems with the experimental values obtained for water and natural
snail mucus. Additionally, if we consider the data given for the natural
logarithm of viscosity as a function of temperature, we see a very similar
trend between the 1MSI-4 system and natural snail mucus. In fact, we know from
Ref. Hanlon  _et al._ (2023a) that for T $\leq$ 290 K, the viscous behaviour
of snail mucus transition from Arrhenius to non-Arrhenius. We can possibly
start to see a similar feature if we consider the single data point at 280 K
($\sim$ 3.6 on x-axis), where the natural logarithm of viscosity increases
substantially at this point (see Fig. 6). Overall, the behaviour observed in
the shear viscosity is consistent with the trends observed in experimental
work, providing strong support that the simulations can accurately describe
the bulk viscoelastic properties of these systems.

##  IV Conclusion

This study is one of very few, if not only, study to systematically
investigate the temperature dependence behavior of viscoelastic properties in
aqueous protein solutions using molecular dynamics simulations, while also
considering variations in protein sizes and concentration. We determined the
temperature dependence of the adiabatic bulk modulus, sound speed, and
viscosity, and found behaviour consistent with that observed in previous
experimental work on protein-water systems Bailey  _et al._ (2020); Comez  _et
al._ (2012); Lupi  _et al._ (2011); Pochylski  _et al._ (2005); Hanlon  _et
al._ (2023a); Trosel  _et al._ (2023); Milewska  _et al._ (2003); Hanlon  _et
al._ (2023b). The general agreement between the simulations results and
experimental observations highlights the potential of the former in predicting
the viscoelastic properties of macromolecular solutions and in providing
complementary molecular-level insight into the behaviour of these systems.

## References

  * Hanlon  _et al._ (2023a) D. F. Hanlon, M. J. Clouter, and G. T. Andrews, “Temperature dependence of the viscoelastic properties of a natural gastropod mucus by brillouin light scattering spectroscopy,” (2023a), arXiv:2306.06807 [cond-mat.mtrl-sci] .
  * Denny (1984) M. Denny, American Zoology 24, 23 (1984).
  * Sharma  _et al._ (2011) V. Sharma, A. Jaishankar, Y.-C. Wang, and G. H. McKinley, Soft Matter 7, 5150 (2011).
  * Dmitriev  _et al._ (2019) A. Dmitriev, V. Vashchenkov, A. Fedoseev, and S. Lushnikov, Journal of Raman Spectroscopy 50, 537 (2019).
  * Hanlon  _et al._ (2023b) D. F. Hanlon, M. J. Clouter, and G. T. Andrews, “Influence of hydration and dehydration on the viscoelastic properties of snail mucus by brillouin spectroscopy,” (2023b), arXiv:2308.00524 [cond-mat.mtrl-sci] .
  * Gautieri  _et al._ (2012) A. Gautieri, S. Vesentini, A. Redaelli, and M. J. Buehler, Matrix Biology 31, 141 (2012).
  * Wei  _et al._ (2010) H. Wei, Y. Fan, and Y. Q. Gao, The Journal of Physical Chemistry B 114, 557 (2010).
  * Silver  _et al._ (2001a) F. H. Silver, J. W. Freeman, I. Horvath, and W. J. Landis, Biomacromolecules 2, 750 (2001a).
  * Silver  _et al._ (2001b) F. H. Silver, D. L. Christiansen, P. B. Snowhill, and Y. Chen, Journal of Applied Polymer Science 79, 134 (2001b).
  * Jia  _et al._ (1996) Z. Jia, C. I. DeLuca, H. Chao, and P. L. Davies, Nature 384, 285 (1996).
  * Durand  _et al._ (2012) E. Durand, E. Derrez, G. Audoly, S. Spinelli, M. Ortiz-Lombardia, D. Raoult, E. Cascales, and C. Cambillau, Journal of Biological Chemistry 287, 38190 (2012).
  * Bujacz (2012) A. Bujacz, Acta Crystallographica Section D: Biological Crystallography 68, 1278 (2012).
  * Hess  _et al._ (2008) B. Hess, C. Kutzner, D. Van Der Spoel, and E. Lindahl, Journal of Chemical Theory and Computation 4, 435 (2008).
  * Lindahl  _et al._ (2001) E. Lindahl, B. Hess, and D. Van Der Spoel, Molecular Modeling Annual 7, 306 (2001).
  * Van Der Spoel  _et al._ (2005) D. Van Der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark, and H. J. Berendsen, Journal of Computational Chemistry 26, 1701 (2005).
  * Abascal and Vega (2005) J. L. Abascal and C. Vega, The Journal of Chemical Physics 123, 234505 (2005).
  * González  _et al._ (2016) M. A. González, C. Valeriani, F. Caupin, and J. L. Abascal, The Journal of Chemical Physics 145 (2016).
  * Yu  _et al._ (2023) Z. Yu, R. Shi, and H. Tanaka, The Journal of Physical Chemistry B 127, 3452 (2023).
  * Jorgensen  _et al._ (1996) W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, Journal of the American Chemical Society 118, 11225 (1996).
  * Zangi (2018) R. Zangi, ACS omega 3, 18089 (2018).
  * Cuendet and van Gunsteren (2007) M. A. Cuendet and W. F. van Gunsteren, The Journal of Chemical Physics 127, 184102 (2007).
  * Humphrey  _et al._ (1996) W. Humphrey, A. Dalke, and K. Schulten, Journal of Molecular Graphics 14, 33 (1996).
  * Ball (2008) P. Ball, Chemical Reviews 108, 74 (2008).
  * Mallamace  _et al._ (2011) F. Mallamace, C. Corsaro, D. Mallamace, P. Baglioni, H. E. Stanley, and S.-H. Chen, The Journal of Physical Chemistry B 115, 14280 (2011).
  * Zhang  _et al._ (2009) Y. Zhang, M. Lagi, D. Liu, F. Mallamace, E. Fratini, P. Baglioni, E. Mamontov, M. Hagen, and S.-H. Chen, The Journal of Chemical Physics 130 (2009).
  * González and Abascal (2010) M. A. González and J. L. Abascal, The Journal of Chemical Physics 132 (2010).
  * Lupi  _et al._ (2011) L. Lupi, L. Comez, C. Masciovecchio, A. Morresi, M. Paolantoni, P. Sassi, F. Scarponi, and D. Fioretto, Journal of Chemical Physics 134, 02B603 (2011).
  * Comez  _et al._ (2012) L. Comez, L. Lupi, M. Paolantoni, F. Picchiò, and D. Fioretto, Journal of Chemical Physics 137 (2012).
  * Pochylski  _et al._ (2005) M. Pochylski, F. Aliotta, Z. Błaszczak, and J. Gapiński, Journal of Physical Chemistry B 109, 4181 (2005).
  * Bailey  _et al._ (2020) M. Bailey, M. Alunni-Cardinali, N. Correa, S. Caponi, T. Holsgrove, H. Barr, N. Stone, C. P. Winlove, D. Fioretto, and F. Palombo, Science Advances 6 (2020).
  * Trosel  _et al._ (2023) Y. Trosel, L. P. Gregory, V. K. Booth, and A. Yethiraj, Biomacromolecules 24, 2469 (2023).
  * Milewska  _et al._ (2003) A. Milewska, J. Szydlowski, and L. P. Rebelo, Journal of Polymer Science Part B: Polymer Physics 41, 1219 (2003).