# Assessing transfer entropy from biochemical data

Takuya Imaizumi, Nobuhisa Umeki1), Ryo Yoshizawa1), Tomoyuki Obuchi2), Yasushi
Sako1), and Yoshiyuki Kabashima3) Department of Mathematical and Computing
Science, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo
152-8550 Japan
1)Cellular Informatics Laboratory, RIKEN Cluster for Pioneering Research, 2-1
Hirosawa, Wako 351-0198, Saitama, Japan
2)Department of Systems Science, Kyoto University, 36-1 Yoshida-Honmachi,
Sakyo-ku, Kyoto 606-8501, Japan
3)The Institute for Physics of Intelligence, The University of Tokyo, 7-3-1
Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

###### Abstract

We address the problem of evaluating the transfer entropy (TE) produced by
biochemical reactions from experimentally measured data. Although these
reactions are generally non-linear and non-stationary processes making it
challenging to achieve accurate modeling, Gaussian approximation can
facilitate the TE assessment only by estimating covariance matrices using
multiple data obtained from simultaneously measured time series representing
the activation levels of biomolecules such as proteins. Nevertheless, the non-
stationary nature of biochemical signals makes it difficult to theoretically
assess the sampling distributions of TE, which are necessary for evaluating
the statistical confidence and significance of the data-driven estimates. We
resolve this difficulty by computationally assessing the sampling
distributions using techniques from computational statistics. The
computational methods are tested by using them in analyzing data generated
from a theoretically tractable time-varying signal model, which leads to the
development of a method to screen only statistically significant estimates.
The usefulness of the developed method is examined by applying it to real
biological data experimentally measured from the ERBB-RAS-MAPK system that
superintends diverse cell fate decisions. A comparison between cells
containing wild-type and mutant proteins exhibits a distinct difference in the
time evolution of TE while apparent difference is hardly found in average
profiles of the raw signals. Such comparison may help in unveiling important
pathways of biochemical reactions.

††preprint: AIP/123-QED

##  I Introduction

Cells and their processes are regulated by interactions with biomolecules via
chemical reactions. Considerable effort has been made to properly understand
the biomolecules and their interactions, which consititutes a tremendous
amount of knowledge [1]. Particularly in the last two decades, chains or
cascades of reactions have gained increasing attention owing to their
potential role in understanding living things as systems [2]. Various
reactions in cells are mathematically modeled using nonlinear differential
equations and Monte Carlo simulations among others. In addition, large
reaction pathways composed of tens or hundreds of components are drawn as
graphs.

These research efforts have made great progress in understanding the mechanism
of cell function control. However, there are still shortcomings. Although
fundamental reactions are modeled by differential equations precisely,
assessing their significance in the reaction cascade is non-trivial. An
increase in the activation levels of enzymes indicates that certain functions
are emerging. Nevertheless, one cannot quantify the significance of the
reactions by only examining solutions of the differential equations as
relevant activation levels can differ from enzyme to enzyme. In addition,
biochemical reactions are dynamical processes, which indicates that the
temporal alteration of signals is an important medium for transferring
relevant information. Therefore, after the reaction pathways are plotted as
graphs, there is a need to clarify the timing and how large the transmitted
information is to properly control cell function.

For compensating these deficiencies, we focus on the transfer entropy (TE) [3,
4]. TE is a universal measure of the causal relationship between two time
series. It is defined by a functional of the joint distributions of two time
series. This makes it possible to quantify causal significance in a unified
manner, regardless of the physical mechanism that generates the time series,
by representing the generation process in probabilistic models. Some empirical
studies support its effectiveness in detecting and characterizing causal
relations in complex systems [5, 6]. TE has been used for analyzing
information processing in various biological organisms from nervous systems
[7, 8, 9, 10, 11, 12, 13, 14, 15] to gene regulatory networks [16, 17, 18, 19]
owing to its universality and effectiveness.

Generally, there are two ways to assess TE: model-based and data-driven
approaches. In the model-based approach, the generation processes of the time
series are modeled using high-dimensional joint distributions or their
equivalent expressions, in which the TE is evaluated analytically or
numerically. Although this approach can potentially enable the exact
assessment of TE for given models, accurately modeling the generation process
of actual systems is difficult, which practically limits its application range
to the analysis of theoretical models [20]. Meanwhile, in the data-driven
approach, the TE is assessed directly from experimentally measured data, which
avoids the difficulty in accurately modeling actual systems. However, the
data-driven approach requires a sufficient number of simultaneously measured
data to substitute the joint distribution of the objective time series as
samples. Hence, technical constraints have restricted its application mostly
to nervous systems, in which large-scale simultaneous measurement is possible
and the required sample sizes are relatively small because of the stationary
nature of the nervous systems [7, 8, 9, 10, 11, 12, 13, 14, 15]. However, the
recent advancement and development of measurement techniques has greatly
raised the ability of simultaneously measuring the activation levels of
biomolecules with high space and time resolutions [21, 22]. This may make it
possible to assess TE of non-stationary biochemical interactions in cells from
measured data.

In such view on the current situation, we explore the possibilities and
limitations for assessing TE from measured data of biochemical reactions.
Generally, the assessment of TE is a computationally demanding task involving
high-dimensional integrals. For overcoming this difficulty, we employ an
approximation by Gaussian models, in which the assessment can be performed
efficiently by estimating the covariance matrix of the time series. However,
it is difficult to theoretically assess the sampling distributions of TE,
which are indispensable for evaluating statistical confidence and significance
of the data-driven estimates, due to the non-stationary nature of biochemical
signals. We overcome this difficulty by developing a method to screen only
statistically significant estimates based on techniques from computational
statistics, the utility of which is tested by the application to a
theoretically tractable time-varying signal model. In addition, we employ the
method for analyzing data experimentally measured from the ERBB-RAS-MAPK
system of actual cells. Distinct difference in the time alteration of TE is
found in the comparison between cells containing wild-type and mutant
proteins, whereas there is no apparent difference in the average profiles of
raw time series. Such comparison may be useful in identifying unknown pathways
of biochemical reactions.

This paper is organized as follows: we review the definition of TE and how it
can be computed under the Gaussian approximation in Section 2. In addition, we
present methods for assessing the statistical confidence and significance of
the inferred results using techniques of computational statistics. In Section
3, the utility of the methods is examined by applying them to data from a
theoretically tractable time-varying signal model. Based on the results
obtained for the theoretically tractable model, we propose a method for
screening only statistically significant estimates of TE. In Section 4, real-
world data from the ERBB-RAS-MAPK system are examined using the method. The
final section presents the summary and discussion of the study.

##  II Assessment of the transfer entropy from measured data

###  II.1 Transfer entropy

Here, we briefly review the definition of TE [3, 4]. We use a conventional
matrix-vector notation in which the bold type denotes column vectors and the
uppercase type represents matrices or random variables. The symbol $\top$
indicates the matrix-vector transpose, such that the (column) vector
$\bm{x}\in\mathbb{R}^{n}$ is represented as
$\bm{x}=(x_{1},\ldots,x_{n})^{\top}$. Given multiple vectors such as
$\bm{x}\in\mathbb{R}^{n}$, $\bm{y}\in\mathbb{R}^{m}$,
$\bm{z}\in\mathbb{R}^{\ell}$ and $\ldots$, their concatenation is expressed as
$\bm{x}\oplus\bm{y}\oplus\bm{z}\oplus\ldots=(x_{1},\ldots,x_{n},y_{1},\ldots,y_{m},z_{1},\ldots,z_{\ell},\ldots)^{\top}$.
The notation ${\mathcal{N}}(\bm{\mu},\Sigma)$ denotes a Gaussian distribution,
where $\bm{\mu}$ is the mean vector and $\Sigma$ is the covariance matrix.

We suppose two interdependent discrete time stochastic processes $X_{t}$ and
$Y_{t}$ $(t=0,1,\ldots,T)$, and denote
$\bm{X}\equiv(X_{T},X_{T-1},\ldots,X_{0})^{\top}$ and
$\bm{Y}\equiv(Y_{T},Y_{T-1},\ldots,Y_{0})^{\top}$. We also use the notation
$\bm{X}_{t}^{(p)}\equiv(X_{t},X_{t-1},\ldots,X_{t-p+1})^{\top}$ and
$\bm{Y}_{t}^{(q)}\equiv(Y_{t},Y_{t-1},\ldots,Y_{t-q+1})^{\top}$. Let us assume
that $\bm{X}$ and $\bm{Y}$ follow a joint distribution $Q(\bm{x},\bm{y})$.
Under this setup, the uncertainty of $X_{t}$ is assessed by (differential)
entropy

$\displaystyle H(X_{t})=-\int dx_{t}Q(x_{t})\ln Q(x_{t}),$ (1)

and that given $\bm{X}_{t-1}^{(p)}$ is quantified by conditional entropy

$\displaystyle H(X_{t}|\bm{X}_{t-1}^{(p)})=-\int
d\bm{x}_{t-1}^{(p)}dx_{t}Q(\bm{x}_{t-1}^{(p)})Q(x_{t}|\bm{x}_{t-1}^{(p)})\ln
Q(x_{t}|\bm{x}_{t-1}^{(p)}),$ (2)

where $Q(x_{t})$, $Q(\bm{x}_{t-1}^{(p)})$, and $Q(x_{t}|\bm{x}_{t-1}^{(p)})$
represent the marginal and conditional distributions with respect to relevant
variables. All of these distributions are reduced from the joint distribution
$Q(\bm{x},\bm{y})$, and similar notations are employed hereafter without
explanations. $H(X_{t}|\bm{X}_{t-1}^{(p)})$ means the remaining uncertainty of
$X_{t}$ on average given $\bm{X}_{t-1}^{(p)}$. Therefore, the amount of the
reduction of the uncertainty

$\displaystyle
I(X_{t};\bm{X}_{t-1}^{(p)})=H(X_{t})-H(X_{t}|\bm{X}_{t-1}^{(p)}),$ (3)

which is referred to as mutual information between $X_{t}$ and
$\bm{X}_{t-1}^{(p)}$ and can be shown to be non-negative, is interpreted as
the amount of information conveyed from the past $p$ states
$\bm{X}_{t-1}^{(p)}=(X_{t-1},...,X_{t-p})^{\top}$ to the current state
$X_{t}$.

Extending the notion of mutual information, transfer entropy (TE) from $Y$ to
$X$ at time $t$ with lags $q$ and $p$ is defined as

$\displaystyle{\cal T}_{Y\to X}^{(q),(p)}(t)$ $\displaystyle=$ $\displaystyle
I(X_{t};\bm{X}_{t-1}^{(p)},\bm{Y}_{t-1}^{(q)})-I(X_{t};\bm{X}_{t-1}^{(p)})$
(4) $\displaystyle=$ $\displaystyle
H(X_{t}|\bm{X}_{t-1}^{(p)})-H(X_{t}|\bm{X}_{t-1}^{(p)},\bm{Y}_{t-1}^{(q)}),$
(5)

and from $X$ to $Y$ as

$\displaystyle{\cal T}_{X\to Y}^{(p),(q)}(t)$ $\displaystyle=$ $\displaystyle
I(Y_{t};\bm{X}_{t-1}^{(p)},\bm{Y}_{t-1}^{(q)})-I(Y_{t};\bm{Y}_{t-1}^{(q)})$
(6) $\displaystyle=$ $\displaystyle
H(Y_{t}|\bm{Y}_{t-1}^{(q)})-H(Y_{t}|\bm{X}_{t-1}^{(p)},\bm{Y}_{t-1}^{(q)}).$
(7)

Equation (5) indicates how much information is increased regarding $X_{t}$ or
how much uncertainty about $X_{t}$ is reduced by providing
$\bm{Y}_{t-1}^{(q)}$ on top of $\bm{X}_{t-1}^{(p)}$, and similarly for (7).
These mean that TE stands for the significance of past states of one variable
in predicting the current state of the other. In addition, the TE is
asymmetric between $X$ and $Y$ unlike mutual information. Therefore, TE is
employed as a useful measure of information flow that quantifies the causal
relationship between two time series.

###  II.2 TE for Gaussian models

TE assessment is computationally demanding even though it is expressed in
compact forms, such as (5) and (7). This is because the computational
complexity grows exponentially with respect to $p+q+1$ when numerically
evaluating (5) and (7), which can be infeasible even for the minimum lags of
$p=q=1$ in practical situations. However, when $Q(\bm{x},\bm{y})$ is a
multivariate Gaussian distribution, the assessment becomes feasible as entropy
and conditional entropy can be analytically evaluated using covariance
matrices for the Gaussian model.

For explaining this more precisely, we introduce the notation $\Sigma(\bm{U})$
to generally express the covariance matrix of the random vector $\bm{U}$. In
addition, we use the notation $\Sigma(\bm{U},\bm{V})$ to denote the cross-
covariance matrix between $\bm{U}$ and $\bm{V}$, which is composed of their
cross-covariances ${\rm cov}(U_{i},V_{\alpha})$. These provide the partial
covariance [4] of $\bm{U}$ given $\bm{V}\oplus\bm{W}\oplus\ldots$ as

$\displaystyle\Sigma(\bm{U}|\bm{V}\oplus\bm{W}\oplus\ldots)=\Sigma(\bm{U})$
(8)
$\displaystyle\phantom{===}-\Sigma(\bm{U},\bm{V}\oplus\bm{W}\oplus\ldots)\Sigma(\bm{V}\oplus\bm{W}\oplus\ldots)^{-1}.\Sigma(\bm{U},\bm{V}\oplus\bm{W}\oplus\ldots)^{\top}.$
(9)

When $Q(\bm{x},\bm{y})$ is given as a multivariate Gaussian distribution, the
properties of the Gaussian random variables yield a formula

$\displaystyle
H(X_{t}|\bm{X}_{t-1}^{(p)})=\frac{1}{2}\ln\left(\Sigma(X_{t}|\bm{X}_{t-1}^{(p)})\right)+\frac{1}{2}\ln(2\pi{\color[rgb]{0,0,0}e}),$
(10)

where
$\Sigma(X_{t}|\bm{X}_{t-1}^{(p)})=\Sigma(X_{t})-\Sigma(X_{t},\bm{X}_{t-1}^{(p)})\Sigma(\bm{X}_{t-1}^{(p)})^{-1}\Sigma(X_{t},\bm{X}_{t-1}^{(p)})^{\top}$
is the variance of $X_{t}$ conditioned by $\bm{X}_{t-1}^{(p)}$.  The
derivation of this formula is presented in Appendix A.  Similarly, we also
obtain

$\displaystyle
H(X_{t}|\bm{X}_{t-1}^{(p)},\bm{Y}_{t-1}^{(q)})=\frac{1}{2}\ln\left(\Sigma(X_{t}|\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})\right)+\frac{1}{2}\ln(2\pi{\color[rgb]{0,0,0}e}).$
(11)

These formulas provide an expression of TE for the Gaussian time series as

$\displaystyle{\cal T}_{Y\to
X}^{(q),(p)}(t)=\frac{1}{2}\ln\left(\frac{\Sigma(X_{t}|\bm{X}_{t-1}^{(p)})}{\Sigma(X_{t}|\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})}\right),$
(12)

and similarly for ${\cal T}_{X\to Y}^{(p),(q)}(t)$.

Several issues are to be noted with this formula. The first is about the
implications of (12). This formula indicates that TE is determined using only
covariances between the two time series irrespectively of their averages. This
implies that the primary media of information transfer are not the average
profiles of the observed signals but their statistical fluctuations, which may
be counterintuitive. Nevertheless, in the framework of the Gaussian
approximation, the covariances are determined by the Hessian of $-\ln
Q(\bm{x},\bm{y})$ around the averages for general joint distributions
$Q(\bm{x},\bm{y})$. This means that the average profile of signals is not a
unique but still a major factor for determining TE. The second is about the
cost for computation. Given the covariance matrix
$\Sigma(X_{t}\oplus\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})$, the
computational cost for assessing (12) increases as $O((p+q)^{3})$ since the
most computationally intensive part is the assessment of the matrix inversion
$\Sigma(\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})^{-1}$. In most cases, this
is computationally feasible as long as $p$ and $q$ are $O(1)$. The third issue
is that, as mentioned in [4], TE is equivalent to the Granger causality [23],
which has been extensively studied since 1970s, for time series generated by
multi-variate autoregressive (MVAR) models. However, in the framework of the
Granger causality, we must introduce many assumptions about how to describe
the time series by MVAR models, which becomes nontrivial, particularly when
handling non-stationary time series. In contrast, the TE framework is “model
agnostic” [4] as (12) can be assessed directly from the covariances, for which
we need few assumptions. Therefore, the formula of (12) would be more user-
friendly when sufficient knowledge about the data generation process is not
available. The final issue is concerning the validity of resorting to Gaussian
models. The appropriateness of modeling time series as Gaussians may be
criticized for specific physical generation processes. However, in most cases,
assessing the TE from the exact formula of (5) and (7) requires significantly
heavy computations for non-Gaussian models; therefore, Gaussian models are
practically unique choices. Further, fortunately, it is known that even in
non-Gaussian cases non-zero values of (12) imply that non-zero values of the
true TE given by (5) [24]. For these reasons, we use the formula of (12) for
assessing TE.

###  II.3 Assessment of TE from data and its statistical significance

Following the above argument, TE can be assessed from samples of time series
${\mathcal{D}}_{M}=\\{\bm{x}_{\mu},\bm{y}_{\mu}\\}_{\mu=1}^{M}$ by
substituting $\Sigma(X_{t}\oplus\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})$
with its estimate from ${\mathcal{D}}_{M}$ in evaluating (12). A natural
estimator of $\Sigma(X_{t}\oplus\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})$
is

$\displaystyle\hat{\Sigma}(X_{t}\\!\oplus\\!\bm{X}_{t-1}^{(p)}\\!\oplus\\!\bm{Y}_{t-1}^{(q)})\\!=\\!\frac{1}{M-1}\sum_{\mu=1}^{M}(\bm{z}_{t,\mu}^{(p),(q)}\\!-\\!\overline{\bm{z}_{t}^{(p),(q)}})(\bm{z}_{t,\mu}^{(p),(q)}\\!-\\!\overline{\bm{z}_{t}^{(p),(q)}})^{\top},$
(13)

where
$\bm{z}_{t,\mu}^{(p),(q)}\equiv(x_{t,\mu},x_{t-1,\mu},\ldots,x_{t-p,\mu},y_{t-1,\mu},\ldots,y_{t-q,\mu})^{\top}$
represents the concatenation of $\mu$-th samples $x_{t,\mu}$,
$\bm{x}_{t-1,\mu}^{(p)}$, and $\bm{y}_{t-1,\mu}^{(q)}$, and
$\overline{\bm{z}_{t}^{(p),(q)}}\equiv
M^{-1}\sum_{\mu=1}^{M}\bm{z}_{t,\mu}^{(p),(q)}$, respectively. This is a
consistent and unbiased estimator of
$\Sigma(X_{t}\oplus\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})$, which
guarantees that (13) converges to the true covariance as $M\to\infty$ and the
average of (13) with respect to the generation of ${\mathcal{D}}_{M}$ accords
with it for finite $M$. The assessment of (12) using (13) also offers a
consistent estimator of ${\cal T}_{Y\to X}^{(q),(p)}(t)$ under Gaussian
assumptions. However, this is generally biased. This is natural because ${\cal
T}_{Y\to X}^{(q),(p)}(t)$ is a non-negative quantity by nature. Therefore,
even if $X_{t}$ and $\bm{Y}_{t-1}^{(q)}$ are statistically independent,
yielding ${\cal T}_{Y\to X}^{(q),(p)}(t)=0$, the statistical fluctuations in
(13) always results in positive values for the estimates of ${\cal T}_{Y\to
X}^{(q),(p)}(t)$. In addition, the convergence rate of the covariance matrix
estimator of (13) is rather slow [25]. This makes it challenging to
analytically evaluate the sampling distributions of the estimates of TE, which
are indispensable for assessing the statistical confidence and significance of
the inferred results, although it is known that the maximum likelihood
estimator of TE will asymptotically have a $\chi^{2}$-distribution under the
null hypothesis in which the true TE vanishes, in the case of stationary time
series [26, 27].

For practically overcoming this difficulty, we employ computational methods
known as bootstrapping [28]. We construct the sampling distributions in the
following ways depending on the purpose:

#### For confidence interval:

Suppose that an estimate $\hat{{\cal T}}_{Y\to X}^{(q),(p)}(t)$ of (12) is
evaluated for a given dataset ${\mathcal{D}}_{M}$ by substituting
$\Sigma(X_{t}\oplus\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})$ with
$\hat{\Sigma}(X_{t}\oplus\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})$ of (13).
We need to evaluate the degree of statistical fluctuations in $\hat{{\cal
T}}_{Y\to X}^{(q),(p)}(t)$, which is inevitable owing to the finiteness of the
sample size $M$. For this purpose, we handle
$\hat{\Sigma}(X_{t}\oplus\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})$ as if it
is the true covariance matrix and generate a new sample set of size $M$ by
independently drawing samples from an identical distribution
${\mathcal{N}}\left(0,\hat{\Sigma}(X_{t}\oplus\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})\right)$.
This provides a surrogate estimate of TE $\hat{{\cal T}}_{Y\to
X}^{*(q),(p)}(t)$ using
$\hat{\Sigma}^{*}(X_{t}\oplus\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})$,
which is the covariance matrix estimated from the new sample set. We repeat
these procedures for $B(\gg 1)$ times, which results in an empirical
distribution of $\hat{{\cal T}}_{Y\to X}^{*(q),(p)}(t)$. We construct the
$100(1-\alpha)$% $(0<\alpha<1)$ confidence interval (CI) by specifying
$50\alpha$ and $100-50\alpha$ percentile points of the empirical distribution
to evaluate the degree of statistical fluctuation of the estimate. The same
procedure is performed for $\hat{{\cal T}}_{X\to Y}^{*(p),(q)}(t)$ as well.

#### For significance threshold:

When estimating $\hat{{\cal T}}_{Y\to X}^{(q),(p)}(t)$, it always takes a non-
negative value even when the true TE ${{\cal T}}_{Y\to X}^{(q),(p)}(t)$
vanishes. We need to evaluate the distribution of $\hat{{\cal T}}_{Y\to
X}^{(q),(p)}(t)$ for the null hypothesis $H_{0}:{{\cal T}}_{Y\to
X}^{(q),(p)}(t)=0$ to distinguish the obtained value from those of chance
levels. To construct the distribution, we define the covariance matrix
$\Sigma_{H_{0}}(X_{t}\oplus\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})$ for
$H_{0}$ by making all elements of the cross-covariance between
$X_{t}\oplus\bm{X}_{t-1}^{(p)}$ and $\bm{Y}_{t-1}^{(q)}$ vanish in
$\hat{\Sigma}(X_{t}\oplus\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})$, keeping
the other elements fixed. Based on this, we generate a new sample set of size
$M$ by independently drawing samples from an identical distribution
${\mathcal{N}}\left(0,\Sigma_{H_{0}}(X_{t}\oplus\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})\right)$.
This provides an estimate of TE $\hat{{\cal T}}_{Y\to X}^{\\#(q),(p)}(t)$
using
$\hat{\Sigma}^{\\#}(X_{t}\oplus\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})$,
which is the covariance matrix estimated from the new sample set. We repeat
these procedures for $B(\gg 1)$ times, which provides an empirical
distribution of $\hat{{\cal T}}_{Y\to X}^{\\#(q),(p)}(t)$. The $100(1-\alpha)$
$(0<\alpha<1)$ percentile point of the empirical distribution is used as the
significance threshold (ST) to distinguish the estimated value from those of
chance levels with the significance level of $100\alpha\%$. The same procedure
is performed for $\hat{{\cal T}}_{X\to Y}^{\\#(p),(q)}(t)$ as well.

These methods are very simple in terms of technical aspects. The necessary
procedure is just repeating the evaluation of the TE from resampled data many
times. The computational burden for repetition can prevent the execution of
these methods in cases where a considerable amount of computation is required
in obtaining a single estimate. However, in the current case, the
computational cost for assessing TE is only $O((p+q)^{3})$ under the Gaussian
approximation, which would not be an obstacle for the execution. Nevertheless,
the validity of using not the true distribution but the empirical distribution
provided by the estimated covariance matrices for assessing the sampling
distributions, following the plug-in principle [28], is debatable. In the next
section, we examine this issue by applying the methods to a time-varying
signal model, in which the theoretical evaluation of TE is tractable.

##  III Testing the developed method using the theoretically tractable model

As a simple but non-trivial example for which the theoretical assessment of TE
is possible, we consider a $d$-dimensional time-varying state space model
[29], which is expressed as

$\displaystyle\bm{u}_{t}=F_{t}\bm{u}_{t-1}+\bm{\mu}_{t}+\bm{\xi}_{t}^{\rm U},$
(14) $\displaystyle\bm{v}_{t}=G\bm{u}_{t}+\bm{\xi}_{t}^{\rm V}$ (15)

where $t=1,\ldots,T$, and $\bm{u}_{t}\in\mathbb{R}^{d(\geq 2)}$ and
$\bm{v}_{t}\in\mathbb{R}^{2}$, respectively. $\bm{\xi}_{t}^{\rm
U}\in\mathbb{R}^{d}$ and $\bm{\xi}_{t}^{\rm V}\in\mathbb{R}^{2}$ are the
Gaussian noise vectors, both of which are independent in time. Equation (14)
describes how the state vector $\bm{U}_{t}$ evolves in time, being subject to
time-varying parameters $F_{t}\in\mathbb{R}^{d\times d}$ and
$\bm{\mu}_{t}\in\mathbb{R}^{d}$. Meanwhile, (15) defines the measurement
process of $\bm{U}_{t}$. We consider that the first and last components of
$\bm{U}_{t}$ are measured, such that the components of the resulting vector
$\bm{V}_{t}$ are regarded as $X_{t}$ and $Y_{t}$ in Section 2 under the setup
of the measurement matrix $G=\left(\begin{array}[]{ccccc}1&0&\ldots&0&0\cr
0&0&\ldots&0&1\end{array}\right)$.

Two features are noted for (14) and (15):

  * •

Given a set of the time-varying parameters $F_{t}$ and $\bm{\mu}_{t}$, the
resulting time series $\bm{u}_{T},\ldots,\bm{u}_{1}$ and
$\bm{v}_{T},\ldots,\bm{v}_{1}$ are provided as linear combinations of
$\bm{\xi}_{T}^{\rm U},\ldots,\bm{\xi}_{1}^{\rm U}$ and $\bm{\xi}_{T}^{\rm
V},\ldots,\bm{\xi}_{1}^{\rm V}$ added to deterministic time series for a given
initial state $\bm{u}_{0}$. This guarantees that the set of time series
$\bm{U}=\\{\bm{U}_{T},\ldots,\bm{U}_{1}\\}$ and
$\bm{V}=\\{\bm{V}_{T},\ldots,\bm{V}_{1}\\}$ follows a joint Gaussian
distribution that varies over time.

  * •

The mean parameters $\bm{\mu}_{T},\ldots,\bm{\mu}_{1}$ do not contribute to
statistical fluctuations of $\bm{U}$. This indicates that covariances of any
components of $\bm{U}$ and $\bm{V}$ are independent of the mean parameters.

These features enable us to evaluate
$\Sigma(X_{t}\oplus\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})$ and
$\Sigma(Y_{t}\oplus\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})$ efficiently
using recursive equations, which results in the theoretical assessment of TE.
More specifically, the statistical independence of noise vector
$\bm{\xi}_{t}^{\rm U}$ regarding time provides a recursive equation for
computing covariance matrix $C_{t}^{\rm
U}=\mathbb{E}[\bm{U}_{t}\bm{U}_{t}^{\top}]-\mathbb{E}[\bm{U}_{t}]\mathbb{E}[\bm{U}_{t}]^{\top}$,
which is expressed as

$\displaystyle C_{t+1}^{\rm U}=F_{t}C_{t}^{\rm U}F_{t}^{\top}+\Delta_{t}^{\rm
U},$ (16)

where $\mathbb{E}[\cdots]$ is the average with respect to noise vectors and
$\Delta_{t}^{\rm U}$ is the covariance matrix of $\bm{\xi}_{t}^{\rm U}$. In
addition, the covariance matrix of $\bm{U}_{t}$ between two different times
are given as

$\displaystyle D_{t+m,t}^{\rm
U}=\left(\prod_{\tau=t+1}^{t+m}F_{\tau}\right)C_{t}^{\rm U},$ (17)

where $m=1,\ldots,T-t$ and $D_{t+m,t}^{\rm
U}=\mathbb{E}[\bm{U}_{t+m}\bm{U}_{t}^{\top}]-\mathbb{E}[\bm{U}_{t+m}]\mathbb{E}[\bm{U}_{t}]^{\top}$.
Subsequently, the covariance matrices of $\bm{V}_{t}$ are computed as

$\displaystyle C_{t}^{\rm V}=GC_{t}^{\rm U}G^{\top}+\Delta_{t}^{\rm V},$ (18)
$\displaystyle D_{t+m,t}^{\rm V}=GD_{t+m,t}^{\rm U}G^{\top},$ (19)

where $C_{t}^{\rm
V}=\mathbb{E}[\bm{V}_{t}\bm{V}_{t}^{\top}]-\mathbb{E}[\bm{V}_{t}]\mathbb{E}[\bm{V}_{t}]^{\top}$,
$D_{t+m,t}^{\rm
V}=\mathbb{E}[\bm{V}_{t+m}\bm{V}_{t}^{\top}]-\mathbb{E}[\bm{V}_{t+m}]\mathbb{E}[\bm{V}_{t}]^{\top}$,
and $\Delta_{t}^{\rm V}$ is the covariance matrix of $\bm{\xi}_{t}^{\rm V}$.
The recursive equations (16) and (17), in conjunction with (18) and (19),
offer all components that constitute $\Sigma(\bm{X}\oplus\bm{Y})$ with
$O(T^{2})$ computational cost. Performing this is feasible as long as $p$ and
$q$ are $O(1)$. Picking up all components of
$\Sigma(X_{t}\oplus\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})$ from
$\Sigma(\bm{X}\oplus\bm{Y})$ and substituting them into (12) provide ${\cal
T}_{Y\to X}^{(q),(p)}(t)$, and similarly for ${\cal T}_{X\to Y}^{(p),(q)}(t)$.

Figure 1:  Theoretically computed ${\cal T}_{{X}\to{Y}}^{(p),(q)}(t)$ for the
system characterized by (25). From top to bottom: The cascade length $d$
varies from $2$ to $5$. From left to right: The lag $p$ of the cause time
series $X_{t}$ changes from $1$ to $4$. See the text for details.  Figure 2:
Computational graph of (25) in the case of $d=4$.

For simulating time series that are simultaneously measured from a cascade of
reactions, we consider a very simple model, which is given by

$\displaystyle F_{t}=\left(\begin{array}[]{ccccc}\alpha&0&\ldots&\ldots&0\cr
F_{21}(t)&\alpha&0&\ldots&0\cr
0&\beta&\ddots&\ddots&\vdots\cr\vdots&\vdots&\ddots&\alpha&0\cr
0&0&\ddots&\beta&\alpha\end{array}\right),$ (25)

where $0<\alpha<1$ is the decay rate of each component and $\beta$ is the
activation caused by the reaction with the component of the previous
subscript. We assume that the first component of $\bm{U}$ triggers the
reaction cascade in a short duration $\gamma>0$, which is taken into account
by the time-varying matrix element

$\displaystyle
F_{21}(t)=\left\\{\begin{array}[]{ll}\delta\exp(-t/\gamma),&t\geq 0\cr
0,&t<0\end{array}\right..$ (28)

Indeed, we need to employ nonlinear equations with respect to $\bm{u}_{t}$ to
precisely describe realistic cascades of chemical reactions. However, the
media for information transfer in the framework of the Gaussian treatment is
the statistical fluctuations of the state vectors around their averages, which
justifies the current linearized description at least as a first
approximation. In addition, the uniform setting of the parameters is not
crucial when examining how the length of the reaction cascade $d$ is reflected
in the dependence of TE on the lag parameters $p$ and $q$, which we will focus
on in the following discussions.

We set $\alpha=0.5$ and $\beta=1$ to ensure that the influence of the trigger
of the first component propagates forward without decay. Meanwhile, we set
$\delta=5$ and $\gamma=1$, $\Delta_{t}^{\rm U}=0.1^{2}\times I_{d\times d}$,
and $\Delta_{t}^{\rm V}=0$ for other parameters, where $I_{n\times n}$
generally denotes $n\times n$ identity matrix. Figure 1 presents the plot of
${\cal T}_{{X}\to{Y}}^{(p),(q)}(t)$ computed from (16)–(19) for various pairs
of lag parameters $p$ and $q$ by varying the length of the reaction cascade
$d$ from $2$ to $5$. The TE of the reverse direction can also be computed, but
it is not presented because the statistical independence between the first
component of $\bm{U}_{t}$ and the subsequent components in the past state
$\bm{U}_{t-m}$ ($m\geq 1$) guarantees that ${\cal T}_{{Y}\to{X}}^{(q),(p)}(t)$
vanishes trivially. All the plots indicate that ${\cal
T}_{{X}\to{Y}}^{(p),(q)}(t)$ does not significantly depends on the lag $q$ of
the effect time series $Y_{t}$. This is presumably because in the current
setup, $\Delta_{t}^{\rm V}=0$ is set to zero for a relatively small $\alpha$.
When $\Delta_{t}^{\rm V}$ is set finite, the dependence of the TE on the lag
$q$ of the effect time series varies in a non-trivial manner depending on the
value of $\alpha$.

On the other hand, ${\cal T}_{{X}\to{Y}}^{(p),(q)}(t)$ monotonically increases
as the lag $p$ of the cause time series $X_{t}$ increases, and almost
saturates at $p=d-1$. This is because the influence of the first component
$U_{1,t}(=X_{t})$ of $\bm{U}_{t}$ reaches the last component $U_{d,t}(=Y_{t})$
late by the lag of $d-1$ due to to the one-dimensional nature of the reaction
cascade, as shown in the computational graph in Figure 2. This graph indicates
that the prediction accuracy of $Y_{t}$ (red node) is maximized when
$U_{d-1,t-1}$ is given in addition to $Y_{t-1}$. Unfortunately, $U_{d-1,t-1}$
is a hidden variable and cannot be directly observed. However, it is
correlated with past states $X_{t-1}(=U_{1,t-1}),X_{t-2}(=U_{1,t-2}),\ldots$
as they share the same ancestors, and the strength of the correlation is
maximized by $X_{t-d+1}$ as it is connected to $U_{d-1,t-1}$ without decaying
through the path of the sequence of yellow nodes in the graph. Therefore,
${\cal T}_{{X}\to{Y}}^{(p),(q)}(t)$ increases as $p$ increases from $1$ to
$d-1$. However, for $p>d-1$, the information from $X_{t-d},\ldots,X_{t-p}$ has
already been considered in
$Y_{t-1}(=U_{d,t-1}),\ldots,Y_{t-p+d}(=U_{d,t-p+d})$, which are denoted as
light green nodes in the graph, and yields little innovative gain on top of
$\bm{Y}_{t-1}^{(q)}$ in predicting $Y_{t}$. Therefore, ${\cal
T}_{{X}\to{Y}}^{(p),(q)}(t)$ hardly increases for $p>d-1$, which is the reason
for the saturation at $p=d-1$.  Such dependence of TE on the lag parameter of
the cause time series does not change qualitatively unless decay rates
$\alpha$ and $\beta$ vary depending on sites significantly around the above-
mentioned values. This, in conjunction with rough knowledge about the time
scale of elemental reactions, may be useful for characterizing the length of
the reaction cascade.

Figure 3:  (a): ${\cal T}_{{X}\to{Y}}^{(p),(q)}(t)$ assessed from 200 samples
for the $d=4$ case presented in Figure 1. (b): That for ${\cal
T}_{{Y}\to{X}}^{(q),(p)}(t)$. In both panels, blue and green full lines stand
for the estimated and true values, respectively. The true values of ${\cal
T}_{{Y}\to{X}}^{(q),(p)}(t)$ are constantly zero in this setting. Light blue
areas and red dotted lines represent 99%CIs and 1%STs. The lower limit of the
99%CI exceeds the true value of the TE at points marked by red crosses, which
may overestimate TE. However, one can screen statistically significant
estimates (marked by red arrows) by accepting only cases in which the lower
limit of the CI is larger than the 1%ST.

Figure 3 shows ${\cal T}_{{X}\to{Y}}^{(p),(q)}(t)$ assessed from $N_{\rm
samp}=200$ samples generated by (14) and (15) for the $d=4$ case presented in
Figure 1 using the methods discussed in Section 2. Estimates of ${\cal
T}_{{Y}\to{X}}^{(q),(p)}(t)$, the true values of which are constantly zero,
are also plotted for reference. For all plots, the upper limit of 99%CIs
surpassed the true value (green) at every point. However, its lower limit is
also larger than the true value at points where the true values are small
(marked by red crosses). In such cases, the CIs do not cover the true values,
leading to inaccurate estimations.

Cross covariances in the estimated matrices
$\hat{\Sigma}(X_{t}\oplus\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})$ and
$\hat{\Sigma}(Y_{t}\oplus\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})$ are
always non-zero even if two time series $\bm{X}$ and $\bm{Y}$ are
statistically independent. This means that TE assessed from samples is always
biased positively even for statistically independent two time series, which
could lead to a risk of giving false positives in judging the finiteness of
TE. One approach to reduce the risk is to estimate the covariance matrices
extremely accurately by collecting a huge number of samples. However, this is
difficult to carry out in practice. Another approach is to take into account
the positive biases in judging the statistical significance. In Figure 3, the
lower limit of the 99%CI is smaller than the 1%ST at all points that are
marked by the red crosses, indicating that the estimates are not statistically
significant in the worst case. In other words, statistically significant
estimates can be screened under given statistical confidence and significance
levels by accepting only cases where the lower limit of the CI is larger than
ST (marked by red arrows). This is the method for assessing TE that we propose
in this study.

The degrees of freedom of covariance matrices
$\hat{\Sigma}(X_{t}\oplus\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})$ and
$\hat{\Sigma}(Y_{t}\oplus\bm{X}_{t-1}^{(p)}\oplus\bm{Y}_{t-1}^{(q)})$ that are
necessary for assessing TE with lag parameters $p$ and $q$ are
$(p+q+2)(p+q+1)/2$. This is as high as $6$ even for the smallest case of
$p=q=1$, and grows up to 28 for $p=q=3$. As $N_{\rm samp}$ must be
sufficiently larger than the degrees of freedom for accurate estimation of the
matrices, the sample size should be at least hundreds even if the proposed
method is employed.

##  IV Application to data from ERBB-RAS-MAPK system

The proposed method for assessing TE was applied to real biological data
obtained from single living cells.

###  IV.1 ERBB-RAS-MAPK system and simultaneous signal measurement by the
fluorescence microscope

The ERBB-RAS-MAPK system is an intracellular signal transduction network
responsible for cell fate decisions [30]. ERBB is a cell surface receptor
protein activated by small proteins, including epidermal growth factor (EGF),
applied to the cell culture medium. The ERBB activation is recognized by the
GRB2/SOS protein complex in the cytoplasm and is associated with the ERBB at
the cytoplasmic side of the cell membrane (Figure 4). Subsequently, SOS
activates RAS protein on the membrane to induce an association of RAF protein
in the cytoplasm with the active RAS. The recruitment of SOS and RAF to the
cell surface can be detected under a fluorescence microscope, reflecting the
activation of ERBB and RAS, respectively [21]. SOS and RAF can be observed
simultaneously in the same single cells using two different colors of
fluorescent tags [22].

Figure 4:  Signal transduction pathway of the ERBB-RAS-MAPK system. In cells
stimulated with EGF, successive translocations of SOS and RAF from the
cytoplasm to the cell membrane occur, which recognizes ERBB and RAS
activation, respectively. The membrane associations of SOS and RAF and RAF-
induced MAPK activation create a negative feedback loop, which is disrupted by
the R1131K mutation of SOS. See the text for details.  Figure 5:  Recruitments
of SOS and RAF to the cell membrane in cells stimulated with EGF. Fluorescence
signals from the basal surface of the cells were selectively observed using a
total internal reflection fluorescence microscope. At time 0, cells were
stimulated with EGF. The upper (SOS) and lower (RAF) images were acquired at
the same field of view using a dual color microscopy. Scale bar: 10 $\mu$m.
Figure 6:  Top panels: 248 samples of raw signals measured from cells with wt
SOS. Bottom panels: Their averages together with one standard deviation. Left
((a), (c)) and right ((b), (d)) panels correspond to SOS and RAF,
respectively.
Figure 7:  Top panels: 282 samples of raw signals measured from cells with
mutant (R1131K) SOS. Bottom panels: Their averages together with one standard
deviation. Left ((a), (c)) and right ((b), (d)) panels correspond to SOS and
RAF, respectively.

###  IV.2 Experiment and measured signals

We stimulated ERBB in HeLa cells expressing the wild-type or mutant (i.e.,
R1131K) SOS with EGF under a microscope at time $t=0$ and measured the changes
in the fluorescence signals from SOS and RAF on the cell membrane (Figure 5).
R1131K is a mutant of SOS, in which the arginine (R) at position 1131 in the
amino acid sequence is replaced by lysine (K). This mutation is observed in
Noonan syndrome, a human genetic disease [31]. The hyperactivation of RAS has
been reported under this mutation, but its molecular mechanism is not entirely
known. One possibility is the defect of a negative feedback regulation caused
by serine phosphorylation around R1131, which is in the GRB2 association site
of SOS [21, 32].

Figures 6 and 7 show 248 and 282 samples of the signals measured from cells
with wild-type (wt) or mutant SOS, respectively. The unit time is 1 minute,
and the signals represent the increase of the fluorescence intensity from the
average levels of the first three points corresponding to $t=-3,-2$, and $-1$
[min]. As shown in these figures, transient increases of the fluorescence
signals of SOS and RAF are observed after EGF application. In addition, the
RAF activation dynamics are delayed and sustained after the initial peak. The
difference in the activation dynamics between cells with wt and mutant SOS is
barely noticed in the average and single-cell trajectories.

Figure 8:  Assessed TE (blue full lines) for cells with wt SOS from samples
shown in Figure 6. (a): From SOS to RAF. (b): From RAF to SOS. Light blue
areas and red dotted lines represent 99%CIs and 1%STs, respectively. The red
arrows indicate statistically significant estimates, which are screened by the
statistical confidence and the significance levels. In response to the initial
peak of TE from SOS to RAF, the significant TE in the reverse direction was
found after 2–3 min (insets). Statistically significant TE is also observed
even in the later stage.  Figure 9:  Assessed TE (blue full lines) for cells
with mutant SOS from samples shown in Figure 7. (a): From SOS to RAF. (b):
From RAF to SOS. Light blue areas and red dotted lines represent 99%CIs and
1%STs, respectively. The red arrows indicate statistically significant
estimates, which are screened by the statistical confidence and the
significance levels. No statistically significant TE was found in the
direction from RAF to SOS, even in the early stages (insets).

###  IV.3 Detecting the difference between cells with wild-type and mutant
biomolecule by TE

Figures 8 and 9 show the assessment of TE values for the samples presented in
Figures 6 and 7 by changing the lags $p$ and $q$ systematically.  A code and
raw data for reproducing these figures are available from a web site[33]. The
sample sizes of 248 and 282 for cells with wt and mutant SOS, respectively,
are considered not to be too small compared with the degrees of freedom
$(p+q+2)(p+q+1)/2$ of the covariance matrices to be estimated within the range
of $1\leq p\leq 3$ and $1\leq q\leq 3$. For all cases, the number of
bootstrapping repetitions $B$ for assessing CI and ST was set to 1000.

Significantly large values of TE from SOS to RAF are observed at the early
stage of cell signaling in all $(p,q)$ combinations for both wt and mutant
SOS, which is expected from the signal transduction cascade (Figure 4).
However, from RAF to SOS, significant TE values are obtained only in cells
with wt SOS (Figure 8). The initial peaks of TE from RAF to SOS delay by
approximately 2–3 minutes from those from SOS to RAF (insets), suggesting
reversed information flow caused by the negative feedback loop from RAF to SOS
via MAPK and/or other proteins downstream of RAF. In addition, the initial
peak of TE from RAF to SOS increases as the lag $q$ of the RAF is set larger,
while that of the reverse direction does not exhibit such a tendency. The
result of the theoretical model in Section 2 implies that this may be because
there is a longer cascade of reactions in the signal transduction pathway from
RAF to SOS than in that from SOS to RAF. The first peak values of TE from RAF
to SOS decrease as $p$ increases, while those from SOS to RAF hardly depend on
$q$. This may be due to the difference of noise levels in the measurement
between the SOS and RAF. The noise levels in our measurement are affected by
the nature of the fluorescence tag bound to the proteins. We used
tetramethylrhodamine and GFP for SOS and RAF, respectively. The former is a
chemical probe that is brighter and more stable than fluorescent proteins like
GFP.

Another striking feature is that statistically significant TE values are
observed in cells with wt SOS at the later stage, while no such behavior is
found in cells with mutant SOS. As RAF activation is sustained in cells with
SOS of both types (Figures 6 and 7), this difference may be a collateral
evidence of the defect of a negative feedback regulation from RAF to SOS in
cells with mutant SOS.

Under the Gaussian approximation, TE can be computed from covariances of two
time series with given time lags. This may invoke a naive question whether
similar results can also be obtained by more conventional covariance based
analyses. For answering such a question, we plot covariance and Pearson’s
correlation coefficient, which is defined by normalizing the covariance with
product of standard deviations, between SOS and RAF changing time lags in
Figures 10 and 11 for cells with wild-type and mutant SOS, respectively. The
plots show that the two time series with the time lags are correlated with the
statistically significant level in the almost entire range of observation time
regardless of whether the SOS type is wild or mutant. Although the profiles of
peaks and the strength of correlations differ slightly, it is difficult to
find distinct qualitative difference between the two cases from the plots.
This implies that TE is a more suitable measure for characterizing the
information flow conveyed between the two time series with a high time
resolution.

Figure 10:  (a): Covariance and (b): Pearson’s correlation coefficient between
SOS and RAF activities with time lags for cells with wild-type SOS (blue full
lines). 99% CIs (light blue areas) and CTs (red dotted lines), which represent
one percentile point from the top for the null hypothesis, are computed by the
bootstrapping method.  Figure 11:  (a): Covariance and (b): Pearson’s
correlation coefficient between SOS and RAF activities with time lags for
cells with mutant SOS (blue full lines). 99% CIs (light blue areas) and CTs
(red dotted lines), which represent one percentile point from the top for the
null hypothesis, are computed by the bootstrapping method.

##  V  Summary and discussion

In summary, we examined the possibilities and limitations of assessing the
transfer entropy (TE) from the measured data of biochemical reactions. We
employed the Gaussian approximation, which enables us to efficiently assess TE
based on covariance matrices estimated from samples of objective time series.
In general, it is necessary to evaluate the sampling distributions to
guarantee the accuracy of the estimated results. However, an analytical
evaluation of the sampling distributions of TE is difficult for non-stationary
time series. We resolved this difficulty by computationally assessing the
sampling distributions using bootstrapping techniques from computational
statistics. The computational methods were tested by the application to a
theoretically tractable model of a stochastic process, which led to the
development of a method for screening only statistically significant estimates
under given levels of statistical confidence and significance. In addition,
this method was applied to assess the dynamics of the information flow in a
real biological reaction network inside living cells. Although the raw signals
measured from cells with the wild-type and a mutant molecule are hardly
distinguished, the method successfully detected the difference between them in
the time course of TE. This implies that the developed method may serve as a
useful tool in studying intracellular reaction networks in which large scale
simultaneous measurement of activities of biomolecules has been made possible
owing to the recent advancement of fluorescence microscope technologies.

It should be noted that the finiteness of TE is not a sufficient sign of the
causal relationship between two time series [34]; non-zero TE can be a
spurious causality detector when observations are performed incompletely due
to the presence of unobserved states [35]. Nevertheless, the TE based analysis
shown in the current paper would still be useful at least for the following
two purposes. One is to quantify the efficiency of information transmission
for known pathways for which unobserved states are absent or their influence
is negligible. As mentioned in Introduction, many pathways of biochemical
reactions have been identified with a high accuracy for these decades.
However, little is known about when and how large information is transmitted
through the pathways. Our methodology could fill the “missing piece” with a
high time resolution although hot debates are still being continued on the
appropriateness of TE as a causality quantifier [36, 37, 38, 39, 40, 41]. The
other is to offer clues for finding unknown pathways. Although the finiteness
of TE is not a sufficient condition of the causal relationship, it still
serves as a necessary condition. Therefore, the assessment of TE would provide
useful guidelines for screening possible candidates of relevant pathways. In
addition, even if pathways are not identified completely, comparison of TE
between healthy and disordered systems might lead to more efficient diagnosis
and treatment of various disorders.

Although we restricted the application domain to biochemical reactions in this
study, the proposed methodology can be utilized to analyze information flow in
general systems of a wider class.  For instance, it may be useful in examining
time varying effective interactions in multi-agent/non-linear dynamical
systems from simulation data [42, 43, 44] and those in nervous/active matter
systems from video imaging data [45, 46], as phenomena observed in such
systems are fairly reproducible and collecting many data on them is relatively
easy. Nevertheless, the necessity of collecting many samples of simultaneously
measured data (i.e., at least hundreds of samples) for accurately estimating
covariances may be a bottleneck in applying it to non-stationary time series
in many other domains.  Meanwhile, in general, appropriate prior knowledge
about objective systems can improve the estimation accuracy significantly. The
reduction of the necessary sample size incorporating Bayesian inference and/or
other machine learning techniques is an important research direction for
further research.

## Acknowledgments

We acknowledge the technical assistance provided by Mutsumi Nakanishi. This
work was partially supported by the MEXT KAKENHI Grant No. 19H05647 (YS), JSPS
KAKENHI Grant Nos. 17H00764 (TO, YK), 18K11463, 19H0182 (TO), 20H00620 (YK),
and JST. CREST Grant No. JPMJCR1912 (YS, YK).

##  Appendix A Derivation of equation (10)

For two random variables $\bm{X}$ and $\bm{Y}$ that follow joint distribution
$Q(\bm{x},\bm{y})$, where $\bm{X}$ and $\bm{Y}$ may be either scalar or any
dimensional vector random variables, conditional entropy of $\bm{X}$ given
$\bm{Y}$ is defined generally as

$\displaystyle H(\bm{X}|\bm{Y})=-\int
d\bm{y}Q(\bm{y})d\bm{x}Q(\bm{x}|\bm{y})\ln Q(\bm{x}|\bm{y}).$ (29)

Let us suppose that $\bm{X}$ and $\bm{Y}$ follows a multivariate Gaussian,
i.e.,

$\displaystyle Q(\bm{x},\bm{y})=\frac{1}{(2\pi)^{\frac{n+m}{2}}({\rm
det}\Sigma(\bm{X}\oplus\bm{Y}))^{1/2}}$ (30)
$\displaystyle\phantom{Q(\bm{x},\bm{y})=}\times\exp\left(-\frac{1}{2}(\bm{x}\oplus\bm{y}-\bm{\mu}_{X}\oplus\bm{\mu}_{Y})^{\top}\Sigma(\bm{X}\oplus\bm{Y})^{-1}(\bm{x}\oplus\bm{y}-\bm{\mu}_{X}\oplus\bm{\mu}_{Y})\right),$
(31)

where $n$ and $m$ are the dimensions of $\bm{X}$ and $\bm{Y}$, and
$\bm{\mu}_{X}$ and $\bm{\mu}_{Y}$ are means of $\bm{X}$ and $\bm{Y}$.
Covariance matrix $\Sigma(\bm{X}\oplus\bm{Y})$ is expressed as

$\displaystyle\Sigma(\bm{X}\oplus\bm{Y})=\left(\begin{array}[]{cc}\Sigma(\bm{X})&\Sigma(\bm{X},\bm{Y})\cr\Sigma(\bm{Y},\bm{X})&\Sigma(\bm{Y})\end{array}\right)$
(34)

using the notation defined in the main text.

Computing the matrix inversion of (34) yields an expression

$\displaystyle\Sigma(\bm{X}\\!\\!\oplus\\!\\!\bm{Y})^{-1}=\\!\\!\left(\\!\\!\begin{array}[]{cc}\Sigma(\bm{X}|\bm{Y})^{-1}&-\Sigma(\bm{X}|\bm{Y})^{-1}\Sigma(\bm{X}\\!,\\!\bm{Y})\Sigma(\bm{Y})^{-1}\cr-\Sigma(\bm{Y})^{-1}\Sigma(\bm{Y}\\!,\\!\bm{X})\Sigma(\bm{X}|\bm{Y})^{-1}&\Sigma(\bm{Y}|\bm{X})^{-1}\end{array}\\!\\!\right)$
(37)

This means that conditional distribution $Q(\bm{x}|\bm{y})$ is expressed as

$\displaystyle Q(\bm{x}|\bm{y})=\frac{1}{(2\pi)^{n/2}{\rm
det}(\Sigma(\bm{X}|\bm{Y}))^{1/2}}$ (38)
$\displaystyle\phantom{Q(\bm{x}|\bm{y})=}\times\exp\left(-\frac{1}{2}(\bm{x}-\bm{\mu}_{X|Y})^{\top}\Sigma(\bm{X}|\bm{Y})^{-1}(\bm{x}-\bm{\mu}_{X|Y})\right),$
(39)

where
$\bm{\mu}_{X|Y}=\bm{\mu}_{X}+\Sigma(\bm{X},\bm{Y})\Sigma(\bm{Y})^{-1}(\bm{y}-\bm{\mu}_{Y})$.
Inserting this expression into (29) offers

$\displaystyle H(\bm{X}|\bm{Y})=\frac{1}{2}\ln{\rm
det}\left(\Sigma(\bm{X}|\bm{Y})\right)+\frac{n}{2}\ln(2\pi e).$ (40)

Finally, substituting $X_{t}$ and $\bm{X}_{t-1}^{(p)}$ with $\bm{X}$ and
$\bm{Y}$, respectively, in (40) provides (10).

## References

  * Alberts  _et al._ [2002] B. Alberts, D. Bray, J. Lewis, M. Raff, K. Roberts, and J. Watson,  _Molecular Biology of the Cell_ , 4th ed. (Garland, 2002).
  * Tavassoly, Goldfarb, and Iyengar [2018] I. Tavassoly, J. Goldfarb, and R. Iyengar, “Systems biology primer: the basic methods and approaches,” Essays in Biochemistry 62, 487–500 (2018), https://portlandpress.com/essaysbiochem/article-pdf/62/4/487/887339/ebc-2018-0003c.pdf .
  * Schreiber [2000] T. Schreiber, “Measuring information transfer,” Phys. Rev. Lett. 85, 461–464 (2000).
  * Barnett, Barrett, and Seth [2009] L. Barnett, A. B. Barrett, and A. K. Seth, “Granger causality and transfer entropy are equivalent for gaussian variables,” Phys. Rev. Lett. 103, 238701 (2009).
  * Bossomaier  _et al._ [2016] T. Bossomaier, L. Barnett, M. Harr, and J. T. Lizier,  _An Introduction to Transfer Entropy: Information Flow in Complex Systems_ , 1st ed. (Springer Publishing Company, Incorporated, 2016).
  * Lungarella  _et al._ [2007] M. Lungarella, K. Ishiguro, Y. Kuniyoshi, and N. Otsu, “Methods for quantifying the causal structure of bivariate time series,” International Journal of Bifurcation and Chaos 17, 903–921 (2007), https://doi.org/10.1142/S0218127407017628 .
  * Wibral  _et al._ [2011] M. Wibral, B. Rahm, M. Rieder, M. Lindner, R. Vicente, and J. Kaiser, “Transfer entropy in magnetoencephalographic data: Quantifying information flow in cortical and cerebellar networks,” Progress in Biophysics and Molecular Biology 105, 80–97 (2011), brainModes: The role of neuronal oscillations in health and disease.
  * Stetter  _et al._ [2012] O. Stetter, D. Battaglia, J. Soriano, and T. Geisel, “Model-free reconstruction of excitatory neuronal connectivity from calcium imaging signals,” PLOS Computational Biology 8, 1–25 (2012).
  * Ursino, Ricci, and Magosso [2020] M. Ursino, G. Ricci, and E. Magosso, “Transfer entropy as a measure of brain connectivity: A critical analysis with the help of neural mass models,” Frontiers in Computational Neuroscience 14, 45 (2020).
  * Honey  _et al._ [2007] C. J. Honey, R. Kötter, M. Breakspear, and O. Sporns, “Network structure of cerebral cortex shapes functional connectivity on multiple time scales,” Proceedings of the National Academy of Sciences 104, 10240–10245 (2007), https://www.pnas.org/content/104/24/10240.full.pdf .
  * Faes, Nollo, and Porta [2011] L. Faes, G. Nollo, and A. Porta, “Information-based detection of nonlinear granger causality in multivariate processes via a nonuniform embedding technique,” Phys. Rev. E 83, 051112 (2011).
  * Stramaglia  _et al._ [2012] S. Stramaglia, G.-R. Wu, M. Pellicoro, and D. Marinazzo, “Expanding the transfer entropy to identify information circuits in complex systems,” Phys. Rev. E 86, 066211 (2012).
  * Sun and Bollt [2014] J. Sun and E. M. Bollt, “Causation entropy identifies indirect influences, dominance of neighbors and anticipatory couplings,” Physica D: Nonlinear Phenomena 267, 49–57 (2014), evolving Dynamical Networks.
  * Novelli  _et al._ [2019] L. Novelli, P. Wollstadt, P. Mediano, M. Wibral, and J. T. Lizier, “Large-scale directed network inference with multivariate transfer entropy and hierarchical statistical testing,” Network Neuroscience 3, 827–847 (2019), https://direct.mit.edu/netn/article-pdf/3/3/827/1092467/netn_a_00092.pdf .
  * Shorten, Spinney, and Lizier [2021] D. P. Shorten, R. E. Spinney, and J. T. Lizier, “Estimating transfer entropy in continuous time between neural spike trains or other event-based data,” PLOS Computational Biology 17, 1–45 (2021).
  * Tung  _et al._ [2007] T. Q. Tung, T. Ryu, K. H. Lee, and D. Lee, “Inferring gene regulatory networks from microarray time series data using transfer entropy,” in  _Twentieth IEEE International Symposium on Computer-Based Medical Systems (CBMS’07)_ (IEEE, 2007) pp. 383–388.
  * Roy  _et al._ [2013] S. Roy, D. Das, D. Choudhury, G. G. Gohain, R. Sharma, and D. K. Bhattacharyya, “Causality inference techniques for in-silico gene regulatory network,” in  _Mining Intelligence and Knowledge Exploration_ , edited by R. Prasath and T. Kathirvalavakumar (Springer International Publishing, Cham, 2013) pp. 432–443.
  * Castro  _et al._ [2019] J. C. Castro, I. Valdés, L. N. Gonzalez-García, G. Danies, S. Cañas, F. V. Winck, C. E. Ñústez, S. Restrepo, and D. M. Riaño-Pachón, “Gene regulatory networks on transfer entropy (grnte): a novel approach to reconstruct gene regulatory interactions applied to a case study for the plant pathogen phytophthora infestans,” Theoretical Biology and Medical Modelling 16, 1–15 (2019).
  * Kim  _et al._ [2020] J. Kim, S. T. Jakobsen, K. N. Natarajan, and K.-J. Won, “TENET: gene network reconstruction using transfer entropy reveals key regulatory factors from single cell transcriptomic data,” Nucleic Acids Research 49, e1–e1 (2020), https://academic.oup.com/nar/article-pdf/49/1/e1/35596141/gkaa1014.pdf .
  * Garcia Michel  _et al._ [2021] L. R. Garcia Michel, C. D. Keirns, B. C. Ahlbrecht, and D. A. Barr, “Calculating transfer entropy from variance-covariance matrices provides insight into allosteric communication in erk2,” Journal of Chemical Theory and Computation 17, 3168–3177 (2021), pMID: 33929855, https://doi.org/10.1021/acs.jctc.1c00004 .
  * Nakamura  _et al._ [2017] Y. Nakamura, N. Umeki, M. Abe, and Y. Sako, “Mutation-specific mechanisms of hyperactivation of noonan syndrome sos molecules detected with single-molecule imaging in living cells,” Scientific Reports 7, 14153 (2017).
  * Yoshizawa  _et al._ [2021] R. Yoshizawa, N. Umeki, A. Yamamoto, M. Okada, M. Murata, and Y. Sako, “p52shc regulates the sustainability of erk activation in a raf-independent manner,” Molecular Biology of the Cell 32, 1838–1848 (2021), pMID: 34260260, https://doi.org/10.1091/mbc.E21-01-0007 .
  * Granger [1969] C. W. J. Granger, “Investigating Causal Relations by Econometric Models and Cross-Spectral Methods,” Econometrica 37, 424–438 (1969).
  * Marinazzo, Pellicoro, and Stramaglia [2008] D. Marinazzo, M. Pellicoro, and S. Stramaglia, “Kernel method for nonlinear granger causality,” Phys. Rev. Lett. 100, 144103 (2008).
  * Ledoit and Wolf [2004] O. Ledoit and M. Wolf, “A well-conditioned estimator for large-dimensional covariance matrices,” Journal of Multivariate Analysis 88, 365–411 (2004).
  * Granger [1963] C. Granger, “Economic processes involving feedback,” Inf. Control. 6, 28–48 (1963).
  * Whittle [1953] P. Whittle, “The analysis of multiple stationary time series,” Journal of the Royal Statistical Society. Series B (Methodological) 15, 125–139 (1953).
  * Efron and Tibshirani [1993] B. Efron and R. J. Tibshirani,  _An Introduction to the Bootstrap_ , Monographs on Statistics and Applied Probability No. 57 (Chapman & Hall/CRC, Boca Raton, Florida, USA, 1993).
  * Durbin and Koopman [2001] J. Durbin and S. J. Koopman,  _Time Series Analysis by State Space Methods_, OUP Catalogue No. 9780198523543 (Oxford University Press, 2001).
  * Lemmon and Schlessinger [2010] M. A. Lemmon and J. Schlessinger, “Cell signaling by receptor tyrosine kinases,” Cell 141, 1117–1134 (2010).
  * Lepri  _et al._ [2011] F. Lepri, A. De Luca, L. Stella, C. Rossi, G. Baldassarre, F. Pantaleoni, V. Cordeddu, B. J. Williams, M. L. Dentici, V. Caputo, S. Venanzi, M. Bonaguro, I. Kavamura, M. F. Faienza, A. Pilotta, F. Stanzial, F. Faravelli, O. Gabrielli, B. Marino, G. Neri, M. C. Silengo, G. B. Ferrero, I. Torrrente, A. Selicorni, L. Mazzanti, M. C. Digilio, G. Zampino, B. Dallapiccola, B. D. Gelb, and M. Tartaglia, “Sos1 mutations in noonan syndrome: molecular spectrum, structural insights on pathogenic effects, and genotype-phenotype correlations,” Human Mutation 32, 760–772 (2011), https://onlinelibrary.wiley.com/doi/pdf/10.1002/humu.21492 .
  * Corbalan-Garcia  _et al._ [1996] S. Corbalan-Garcia, S. S. Yang, K. R. Degenhardt, and D. Bar-Sagi, “Identification of the mitogen-activated protein kinase phosphorylation sites on human sos1 that regulate interaction with grb2,” Molecular and Cellular Biology 16, 5674–5682 (1996).
  * [33] https://github.com/kabashiy/transfer_entropy/.
  * Baldovin, Cecconi, and Vulpiani [2020] M. Baldovin, F. Cecconi, and A. Vulpiani, “Understanding causation via correlations and linear response theory,” Phys. Rev. Research 2, 043436 (2020).
  * Smirnov [2013] D. A. Smirnov, “Spurious causalities with transfer entropy,” Phys. Rev. E 87, 042917 (2013).
  * Smirnov [2014] D. A. Smirnov, “Quantification of causal couplings via dynamical effects: A unifying perspective,” Phys. Rev. E 90, 062921 (2014).
  * Stokes and Purdon [2017] P. A. Stokes and P. L. Purdon, “A study of problems encountered in granger causality analysis from a neuroscience perspective,” Proceedings of the National Academy of Sciences 114, E7063–E7072 (2017), https://www.pnas.org/content/114/34/E7063.full.pdf .
  * Dhamala  _et al._ [2018] M. Dhamala, H. Liang, S. L. Bressler, and M. Ding, “Granger-geweke causality: Estimation and interpretation,” NeuroImage 175, 460–463 (2018).
  * Barnett, Barrett, and Seth [2018] L. Barnett, A. B. Barrett, and A. K. Seth, “Solved problems for granger causality in neuroscience: A response to stokes and purdon,” NeuroImage 178, 744–748 (2018).
  * Smirnov [2020] D. A. Smirnov, “Transfer entropies within dynamical effects framework,” Phys. Rev. E 102, 062139 (2020).
  * Smirnov [2019] D. A. Smirnov, “Spectral causalities within dynamical effects framework,” EPL (Europhysics Letters) 128, 20006 (2019).
  * Vicsek  _et al._ [1995] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, “Novel type of phase transition in a system of self-driven particles,” Phys. Rev. Lett. 75, 1226–1229 (1995).
  * Andrzejak, Ledberg, and Deco [2006] R. G. Andrzejak, A. Ledberg, and G. Deco, “Detecting event-related time-dependent directional couplings,” New Journal of Physics 8, 6–6 (2006).
  * Cavagna and Giardina [2014] A. Cavagna and I. Giardina, “Bird flocks as condensed matter,” Annual Review of Condensed Matter Physics 5, 183–207 (2014), https://doi.org/10.1146/annurev-conmatphys-031113-133834 .
  * Ota  _et al._ [2021] K. Ota, Y. Oisi, T. Suzuki, M. Ikeda, Y. Ito, T. Ito, H. Uwamori, K. Kobayashi, M. Kobayashi, M. Odagawa, C. Matsubara, Y. Kuroiwa, M. Horikoshi, J. Matsushita, H. Hioki, M. Ohkura, J. Nakai, M. Oizumi, A. Miyawaki, T. Aonishi, T. Ode, and M. Murayama, “Fast, cell-resolution, contiguous-wide two-photon imaging to reveal functional network architectures across multi-modal cortical areas,” Neuron 109, 1810–1824.e9 (2021).
  * Iwasawa, Nishiguchi, and Sano [2021] J. Iwasawa, D. Nishiguchi, and M. Sano, “Algebraic correlations and anomalous fluctuations in ordered flocks of janus particles fueled by an ac electric field,” Phys. Rev. Research 3, 043104 (2021).