# How particular is the physics of the free energy principle?

Miguel Aguilera<EMAIL_ADDRESS>Beren Millidge  Alexander Tschantz
Christopher L. Buckley

###### Abstract

The free energy principle (FEP) states that any dynamical system can be
interpreted as performing Bayesian inference upon its surrounding environment.
In this work, we examine in depth the assumptions required to derive the FEP
in the simplest possible set of systems – weakly-coupled non-equilibrium
linear stochastic systems. Specifically, we explore (i) how general the
requirements imposed on the statistical structure of a system are and (ii) how
informative the FEP is about the behaviour of such systems. We discover that
two requirements of the FEP – the Markov blanket condition (i.e. a statistical
boundary precluding direct coupling between internal and external states) and
stringent restrictions on its solenoidal flows (i.e. tendencies driving a
system out of equilibrium) – are only valid for a very narrow space of
parameters. Suitable systems require an absence of perception-action
asymmetries that is highly unusual for living systems interacting with an
environment. More importantly, we observe that a mathematically central step
in the argument, connecting the behaviour of a system to variational
inference, relies on an implicit equivalence between the _dynamics of the
average_ states of a system with the _average of the dynamics_ of those
states. This equivalence does not hold in general even for linear systems,
since it requires an effective decoupling from the system’s history of
interactions. These observations are critical for evaluating the generality
and applicability of the FEP and indicate the existence of significant
problems of the theory in its current form. These issues make the FEP, as it
stands, not straightforwardly applicable to the simple linear systems studied
here and suggest that more development is needed before the theory could be
applied to the kind of complex systems that describe living and cognitive
processes.

###### keywords:

Free energy principle , Markov blanket , Non-equilibrium systems , Bayesian
inference , Linear stochastic systems

††journal: Physics of Life Reviews

[inst1]organization=School of Engineering and Informatics, University of
Sussex,city=Falmer, Brighton, postcode=BN1 9QJ, country=United Kingdom

[inst2]organization=MRC Brain Network Dynamics Unit, University of
Oxford,city=Oxford, postcode=OX1 3TH, country=United Kingdom

[inst3]organization=Sackler Center for Consciousness Science, University of
Sussex ,city=Falmer, Brighton, postcode=BN1 9QJ, country=United Kingdom

##  1 Introduction

During the last decade, the ‘free energy principle’ (FEP) has become an
influential framework which aims to provide a grand theory promoting a
Bayesian interpretation of living systems [1, 2, 3]. The FEP states that any
self-organizing system (i.e. any dynamical system, and therefore any living or
cognitive entity) equipped with a Markov blanket – a statistical separation
between internal and external states – can be interpreted as performing
Bayesian inference upon the surrounding environment, such that its internal
states come to encode probabilistic beliefs about the external environment [3,
4].

The core claim of the FEP is exceptionally ambitious. It implies that the
dynamics of any pair of coupled systems, under specific conditions about the
interaction of internal and external states, can be described as one system
trying to statistically infer the states of the second system (cf an agent and
its environment). This claim licenses an interpretation of the agent as
performing a basic kind of Bayesian inference and encoding beliefs about the
surrounding environment [5, 6]. Such an equivalence could have a far-reaching
influence on the study of living systems. For example, it could enable
approximate calculations of the dynamics of complex systems in terms of a more
tractable description of the dynamics of their sufficient statistics.
Furthermore, the FEP has been defended based on singular claims about its
explanatory power, suggesting that it reveals novel insights among fundamental
psychological concepts such as memory, attention, value, reinforcement, and
salience [7] and unifies different aspects of motor behaviour and perception,
from retinal stabilization to goal-seeking [8]. In addition, it has been
proposed that the FEP provides a basis for integrating several general brain
theories, including the Bayesian brain hypothesis, neural Darwinism, Hebbian
cell assembly theory, and optimal control and decision theory [9].

The FEP has also inspired theories such as predictive coding [10, 11, 12, 13]
and active inference [14, 15, 16], which offer explanations and models of
brain function and dysfunction [17] through the lens of Bayesian inference,
and have become widely influential in theoretical neuroscience and beyond. For
instance, predictive coding has been proposed to be a biologically plausible
model of cortical function [11, 12, 18], and has been applied to explain
binocular rivalry [19] or attention [20]. Active inference, on the other hand,
has seen substantial use in modelling the behaviour of human or animal
subjects in various paradigms [14, 15], as well as understanding rational
decision-making and behavioural control [21, 22, 23, 24]. Moreover, biological
theories proposing that neurons [25], synapses [26], bacteria [27] or plants
[28] could be explicitly performing Bayesian (variational) inference by
minimizing free energy gradients have been proposed and justified using
explicit appeals to the FEP. While such theories do not entirely depend upon
the validity of the mathematical core of the FEP reviewed here, they
nevertheless derive a substantial amount of their intellectual and rhetorical
support from it. Thus, the fundamental validity of the mathematical framework
of the FEP is of great importance to this large and rapidly increasing
modelling literature.

In a vast body of work that spans over the course of about 15 years, the FEP
has detailed the mathematical steps required to derive its central claims. In
the first phase, an intuitive and heuristic idea of an imperative to minimize
variational free energy was developed [29, 7] based on the need for any
recognizable ‘system’ to maintain itself in a low entropy configuration over
time against dissipative forces trying to push it towards a high entropy
equilibrium state. Later, this heuristic argument and intuition was more
formally related to concepts from stochastic thermodynamics [2, 1],
specifically by bringing in the Markov blanket condition and expressing the
dynamics of the system in terms of a gradient descent on a variational free
energy corresponding to a generative model of the environment of the system.
Finally, the mathematical formulation and some of its arguments have been
recently refined in a new series of publications [3, 4, 30], a process that is
still ongoing.

Figure 1: Growth in publications pertaining the free energy principle.
Publications per year searching for ‘free energy principle’. Source: Google
Scholar. The number of papers has approximately doubled from 2013 (link) to
2017 (link) and from 2017 to 2020 (link), as indicated by the dashed lines.

Despite the extensive literature on the FEP, there are few concrete examples
that apply all the required steps to a specific, well-studied system.
Similarly, it has rarely been explored whether the assumptions (which are
sometimes left unstated) required for deriving the principle hold under the
dynamics we expect from cognitive and living systems. In light of the
increasing rate of publications concerning the FEP (the number of published
papers on the topic doubles every few years, Fig. 1), we believe it is
imperative to ground and test the foundations of the theory in concrete models
to assess the generality and validity of its claims. To this end, we explore a
class of systems defined by stochastic linear differential equations, under a
weak-coupling assumption. Such systems are the simplest possible example that
can display the dynamics required for the FEP, as well as capture non-
equilibrium properties of systems engaged in perception-action cycles (as we
expect from living systems). Moreover, the absence of nonlinear interactions
in this class of systems allows for precise analytic calculations, which
offers an interesting test-bed to examine, in detail, the connection between
non-trivial dynamics and the statistical properties of coupled systems. In
general, due to the special independence properties imposed by the theory, if
the assumptions and steps of the FEP do not hold in such simple systems, we
consider it unlikely that they hold in more complex nonlinear systems where
the dynamics are expected to be more deeply intertwined. Moreover, we observe
that stronger couplings result in higher-order interactions that make
difficult the separation between internal and external states required by the
FEP. Because of this, we expect that the introduction of non-linearities will
in general have a similar effect.

### How general is the free energy principle

We first inspect the generality of the assumptions about the statistical
structure of a system required by the principle. A crucial step to derive the
FEP is to establish a relation between the average flow of change of a system
that interacts with an environment and the gradient of a variational free
energy of a model of this environment. This step relies on specific
assumptions about how perception and action mediate the interaction of the
internal and external states of a system. We aim to explore how general these
assumptions are and whether they can be expected in the dynamical systems
models of living systems.

Perception-action interface. The FEP partitions the states of a system into
external, sensory, active and internal states. Then, the theory assumes that
perception-action cycles involve causal dependencies such that internal and
external states are only mutually influenced through the effect of active and
sensory states [2]. We will refer to this idea as a perception-action
interface (see A). One example of this interface is a cell membrane around a
cell [2], although the formal definition of a perception-action interface does
not require the interface to be a physical boundary. Thus, a perception-action
interface could describe any set of variables mediating between system and
environment as, for example, a combination of retinal activity and oculomotor
states mediating between neural activity in the visual cortices and the
location of an object in the environment [4].

Then, the FEP requires the system to be endowed with a particular statistical
structure with two special properties: a Markov blanket (i.e. conditional
independence between internal and external states) and the absence of
solenoidal couplings.

Markov blanket. The FEP prescribes that variables in a perception-action
interface constitute a Markov blanket [2, 30]. A Markov blanket (see A) is
defined as a set of states (the ‘blanket’) that separates two other sets in a
statistical sense (i.e. they are conditionally independent, given the
blanket). The term was initially introduced in the context of Bayesian
networks or graphs [31], and it is also known as the general Markov condition
[32]. In the FEP, Markov blankets are used for identifying a set of variables
that separate the internal and external states of a system (see [33] for a
detailed study on the specific use of the concept of Markov blanket in the
FEP). Here, we note that Markov blankets can be easily identified in models
defined by directed acyclic Bayesian networks (Fig. 2.A). In these systems, a
sufficient condition for a Markov blanket of variable $\mathbf{x}$ (e.g., an
internal state) is that it contains the parent nodes of $\mathbf{x}$, the
children nodes of $\mathbf{x}$ and the other parents of each children node (in
this case, the minimum Markov blanket is composed by the grey nodes,
$\mathbf{s},\mathbf{a}$, in Fig. 2.A). This is defined as the local Markov
condition [32] (which implies the general Markov condition in the directed
acyclic graphs of Bayesian networks). However, in [2] the FEP suggests that a
Markov blanket arises naturally from the perception-action interface depicted
in Fig. 2.B (although recent works restrict this to the case of an absence of
particular solenoidal flows [30]). However, such a cyclic structure can
generate couplings that propagate beyond causal interactions. In this case,
the local Markov condition does not imply a Markov blanket, as marginalization
over local blanket variables generates new couplings (dashed arrows in Fig.
2.B). This important issue (identified by [34]) contradicts the intuition that
perception-action states always constitute Markov blankets, which we will
explore in following sections.

Figure 2: Markov blankets and causal cyclic models. (A) An example of a Markov
blanket (grey nodes $\mathbf{s},\mathbf{a}$, e.g. perception and action) for a
variable $\mathbf{x}$ (e.g. an internal state) in a directed acyclic graph
representing a Bayesian network. (B) An example of a cyclic causal model where
the Markov blanket is not directly identifiable from the graph structure. In
this case, nodes that meet the local Markov condition (grey nodes) do not
guarantee the conditional independence required for a Markov blanket, as new
couplings might emerge (dashed arrows between $\mathbf{x},\mathbf{y}$). In
each case, the graph at the right represents the equivalent correlations of
the system probability distribution for fixed $\mathbf{s},\mathbf{a}$.

Solenoidal couplings. The second required property is that solenoidal
couplings between internal, external, sensor and motor states are absent. The
idea of solenoidal flows (see A) arises from the separation of the flow of a
dynamical system into two components. The first component is a dissipative
(curl-free) flow that counters the dispersion of the density caused by random
fluctuations in the system. The second component, the solenoidal flow, is
defined as a conservative (divergent-free) flow capturing dissipative
tendencies in the system, driving it away from equilibrium [3, p.11]. This
nonequilibrium nature is a fundamental aspect of living entities. Examples of
this are asymmetric organism-environment interactions [35] or the oscillatory
behaviour underpinning most biorhythms and neural dynamics [36]. The FEP
assumes that solenoidal couplings between internal and external states are
absent, meaning that these flows will not penetrate through the perception-
action interface.

Average flows and the variational free energy. The confluence of these
properties entails an important result: due to this particular statistical
structure, the average flow of a system can be described in terms of a
variational free energy gradient. The average flow (also called the marginal
flow, see A) describes the average rate of change of the system conditioned on
the blanket state. In turn, the variational free energy (see A) represents an
upper bound of the surprise of observed states, according to an internal model
of the environment.

In the first section of this manuscript, we study how likely these conditions
are for the type of dynamical interactions we expect to find in living system.
By considering the simplest case of non-equilibrium stochastic dynamics, we
show that Markov blankets and absence of solenoidal flows only emerge for very
particular perception-action interfaces, forcing symmetries in agent-
environment interactions that are not expected in living beings.

The core issue is that these three requirements – a perception-action
partition where sensor and active states mediating the internal and the
external states, the existence of a Markov blanket, and decoupled solenoidal
flows – are, in principle, independent conditions. A perception-action
interface does not necessarily guarantee the required conditional independence
relationships and, in fact, generally does not since statistical correlations
can propagate beyond this interface over time due to the intrinsic
fluctuations and reentrant connections in the system. Conversely, the
conditional independence relationship decreed by the Markov blanket does not
imply that there is a lack of dynamical coupling between internal and external
states [34]. In practice, we discover that the kinds of systems that can
fulfil both the Markov blanket condition and the block-diagonal solenoidal
coupling condition are extremely specialized and generally do not possess the
kind of sparsity and asymmetry of dynamical couplings that we expect from a
perception-action interface. These asymmetries are present at different levels
of living systems, leading to qualitative differences between the inside and
outside of a system, shaping system-environment interactions and its related
flows of energy and matter [35, 37, 38, 39].

In other words, although the Markov blanket assumption is typically maintained
in systems with very weak couplings, we observe that there are direct
interactions between external and internal states in the system’s dynamics.
This subtle distinction between the conditional independence relationships of
the mechanisms of a system (perception-action interface) and its statistical
couplings (Markov blanket) – which is analogous to the distinction between
anatomical and functional connectivity in neuroscience – has perhaps been
underappreciated in the FEP literature. This ambiguity leads to the claim that
the requirements of the FEP naturally describe systems with a causal boundary
between the external and internal states [2], when this is not necessarily the
case.

### How informative is the free energy principle

Once the relation between a free energy gradient and the average flow of a
system has been established, we explore how informative this relationship is
about the behaviour of an organism or the evolution of a dynamical system.

Conditional synchronisation manifold. To justify the relation between a
gradient of a free energy functional and the behaviour of a system, the FEP
assumes the existence of a conditional synchronisation manifold (see A). This
manifold is defined as a mapping that, given a blanket state, connects the
most likely internal and external states. The FEP proposes that its existence
allows us to characterise the relationship between (maximum a posteriori)
internal and external states in terms of internal states ‘sensing’ or
‘tracking’ external states through the Markov blanket [3].

Next, the FEP links the evolution of the most likely external states with the
average flow of the system, conditioned on the blanket states. Since the step
described in the previous section connects the free energy gradient with the
average flow, this assumption implies that the evolution of the most likely
external states is driven by the gradient of a variational free energy.
Moreover, given the conditional synchronisation manifold, the most likely
internal states are also driven by this free energy gradient. The FEP suggests
this result leads to the appearance that any system with the properties
described above behaves _as if_ internal states were performing variational
inference to predict external states.

Average flows and the rate of change of the average. The claim that systems
under the required properties behave as if performing inference relies on a
venturesome assumption. It implicitly attributes the rate of change of the
average (expected) state of a system– i.e. the expected change in the internal
state of an agent at a particular moment – can be roughly described by the
steady-state average flow of the system conditioned on a blanket state – i.e.
the average rate of that state during many trajectories (see A). This
assumption is driven by the intuition that if the average flow points in a
direction that minimizes variational free energy, internal states will behave
as if they are trying, on average, to minimize this free energy. As well, this
can be read as an implicit assumption that these two quantities are
approximately the same. However, as we will show, this intuition is incorrect
since the average flow conditioned on blanket states disconnects the rate of
change in the system from its previous trajectory, which is crucial to predict
the system’s behaviour. In practice, this assumption implies decoupling the
actions of an agent from its history of previous states. We will show that, in
the class of linear system explored in this work, this results in the free
energy gradients being uninformative about the behaviour of an agent or its
specific trajectories.

We should note that this claim has been relaxed in more recent work [30],
proposing (instead of an equivalence between rates of change of expected
states and average rates) that an interpretation in terms of Bayesian
inference emerges only in expectation. We will see however that this still
presents important practical and conceptual problems. To draw an analogy that
portrays this issue, we could propose that the actions of a population of
organisms in a particular evolutionary context maximize, on average, a fitness
function – e.g. the number of genes the population passes on to the next
generation. That is, however, a largely uninformative statement for describing
the behaviour of an individual organism, which depends on that organism’s
specific history. Moreover, the behaviour of an organism that systematically
inferred what to do next in terms of (evolutionary) fitness maximization would
likely be entirely different from the realized behaviour of any living
organism. This is not only a point of philosophical nuance but, as we will
show, translates into assumptions and important equations required for
deriving the FEP.

### Overview of the mathematical review

In this article, we present a technical and conceptual critique of the FEP. As
the theory is in continuous development, its mathematical details have been
described in different forms and notations in the literature. For our study we
try to remain as close as possible to the verbal and formal descriptions
described in the most recent publications [3, 4, 30]. However, in at least two
instances (see Assumption 1 and Assumptions 3∗ and 3∗∗), we identify
conflicting interpretations of the theory. In these cases, we attempt to
derive our argument with as much mathematical coherence as possible while also
presenting results that address the different identified possibilities.

The rest of the article is organized as follows. First, we present a summary
of the FEP, including a list of conditions and assumptions required to derive
it. Next, we survey the steps prescribed to derive the FEP for a linear
stochastic system under a weak coupling assumption. We then evaluate in which
cases we expect the requirements of the FEP to hold. Finally, we evaluate the
implications of our study for the theory and its applicability to the type of
processes that living systems are expected to manifest.

Figure 3: Conditions for the free energy principle. List of conditions and
assumptions for deriving the FEP. The conditions for the FEP are: the flow of
the system is constrained to define a perception-action interface (Condition
1, illustrated in the figure for an Euler integration step), the system has a
global attractor described by a state-independent SDE decomposition (Condition
1), and from this configuration a Markov blanket emerges (Condition 3). The
assumptions necessary for deriving the FEP are: the global attractor presents
an uncoupling of solenoidal flows (Assumption 1), there is an invertible
mapping between the most likely internal and external states (Assumption 2)
and the evolution of the most likely external dynamics is described by the
marginal flow (Assumption 3∗).

##  2 Summary of the theory

Here we present a succinct description of the theory, based on the most recent
publications [3, 4, 30], although some steps apply to previous versions as
well. Generally, the FEP assumes a random dynamical system described in terms
of a Langevin stochastic differential equation:

$\displaystyle\frac{\mathrm{d}\mathbf{z}_{t}}{\mathrm{d}t}=f(\mathbf{z}_{t})+\boldsymbol{\omega}_{t},$
(1)

where $\mathbf{z}_{t}=\\{z_{i,t}\\},i=1\ldots N$ is a vector, $f$ is an
arbitrary but differentiable function and
$\boldsymbol{\omega}_{t}=\\{\omega_{i,t}\\}$ is a Gaussian white noise with
covariance $2\boldsymbol{\Gamma}$, which is a diagonal matrix. Throughout this
article, bold symbols represent vectors and matrices.

The FEP further assumes that the system can be decomposed into external,
sensory, active and internal states,
$\mathbf{z}=\\{\mathbf{y},\mathbf{s},\mathbf{a},\mathbf{x}\\}$, configured as
a perception-action loop reflecting an interface mediating between
‘autonomous’ states (active and internal states $\\{\mathbf{a},\mathbf{x}\\}$)
and ‘non-autonomous’ states (external and sensory states
$\\{\mathbf{y},\mathbf{s}\\}$). This leads to

######  Condition 1.

The flow function $f$ decouples autonomous and non-autonomous states according
to the following perception-action interface:

$\displaystyle
f(\mathbf{z}_{t})=\begin{bmatrix}f_{y}(\mathbf{y}_{t},\mathbf{s}_{t},\mathbf{a}_{t})\\\
f_{s}(\mathbf{y}_{t},\mathbf{s}_{t},\mathbf{a}_{t})\\\
f_{a}(\mathbf{s}_{t},\mathbf{a}_{t},\mathbf{x}_{t})\\\
f_{x}(\mathbf{s}_{t},\mathbf{a}_{t},\mathbf{x}_{t})\end{bmatrix}.$ (2)

Under the presence of random fluctuations, some systems will converge toward a
stable global attractor reflecting the steady-state dynamics of the system. In
systems out of thermodynamic equilibrium, this global attractor will describe
a non-equilibrium steady-state (NESS), characterized by a continuous energy
flux between the system and its environment. The next condition of the FEP is
that this global attractor exists and can be described by a stochastic
differential equation (SDE) decomposition [40, 41] (often referred in the FEP
literature as a ‘Helmholtz’ decomposition) describing the flow function as a
linear function of the logarithmic steady state distribution of the system:

######  Condition 2.

The FEP assumes that the system will reach a non-equilibrium steady state
described by the probability density function $p(\mathbf{z}_{t})$, which can
be described using a SDE decomposition that separates the flow into
dissipative
$(-\boldsymbol{\Gamma}\nabla_{\mathbf{z}}\textfrak{J}(\mathbf{z}_{t}))$ and
solenoidal $(\mathbf{Q}\nabla_{\mathbf{z}}\textfrak{J}(\mathbf{z}_{t}))$
components

$\displaystyle f(\mathbf{z}_{t})=$
$\displaystyle(\mathbf{Q}-\boldsymbol{\Gamma})\nabla_{\mathbf{z}}\textfrak{J}(\mathbf{z}_{t}),$
(3) $\displaystyle\textfrak{J}(\mathbf{z}_{t})=$ $\displaystyle-\log
p(\mathbf{z}_{t}),$ (4)

where $\mathbf{Q}$ is an antisymmetric matrix – i.e. equal to its negative
transpose, $\mathbf{Q}=-\mathbf{Q}^{\intercal}$. This condition requires
$\boldsymbol{\Gamma}$ and $\mathbf{Q}$ to be constant matrices – i.e. state-
independent, as it is the case in linear systems111Recent works have relaxed
this assumption, considering the case of state-dependent $\boldsymbol{\Gamma}$
and $\mathbf{Q}$ matrices [42, 43]. However, these changes do not materially
affect our critique (which uses systems where these matrices are independent
of the state of the system) and we focus on the state-independent case here..

The third condition is that the perception-action interface induces a Markov
blanket into the NESS probability distribution. In a Markov blanket, internal
states are independent when conditioned on the blanket states, composed of
sensory and active states, $\mathbf{b}=\\{\mathbf{s},\mathbf{a}\\}$:

######  Condition 3.

The steady state distribution is described in terms of a Markov blanket, where
internal/external states are independent when conditioned on its blanket
states.

$\displaystyle
p(\mathbf{y}_{t},\mathbf{x}_{t}|\mathbf{b}_{t})=p(\mathbf{y}_{t}|\mathbf{b}_{t})p(\boldsymbol{x}_{t}|\mathbf{b}_{t}).$
(5)

It is important to note (see [34]), that Conditions 1 and Condition 3 are
independent and one does not entail the other. This is an important point, as
some presentations of the FEP assume that a perception-action interface
directly involves a causal barrier (Markov blanket) between system and
environment, and this is not always the case.

###  2.1 First move: Capturing Bayesian inference with an average flow

The FEP starts by describing the average flows of the external states of a
system as following a gradient minimizing a variational free energy. This
connects these flows with notions from Bayesian inference.

The principle starts from a description of the ‘surprise’ of the observed
blanket states, under the steady-state distribution defined by the (random)
dynamics of the system, which we denote with $\textfrak{J}(\mathbf{b})$ being
the negative log-probability

$\displaystyle\textfrak{J}(\mathbf{b})=-\log p(\mathbf{b}).$ (6)

This implies that highly unlikely states will have a large surprise value and
vice versa.

However, a system cannot access this surprise value without complete knowledge
of its environment. Thus, Bayesian inference prescribes to use a lower bound
of this surprise described by the variational free energy
$F(\boldsymbol{\theta},\mathbf{b})$:

$\displaystyle-\log p(\mathbf{b})<F(\boldsymbol{\theta},\mathbf{b})=$
$\displaystyle-\log
p(\mathbf{b})+D_{KL}(q(\mathbf{y}|\boldsymbol{\theta})||p(\mathbf{y}|\mathbf{b}))$
(7) $\displaystyle
D_{KL}(q(\mathbf{y}|\boldsymbol{\theta})||p(\mathbf{y}|\mathbf{b}))=$
$\displaystyle\int
d\mathbf{y}q(\mathbf{y}|\boldsymbol{\theta})\log\frac{q(\mathbf{y}|\boldsymbol{\theta})}{p(\mathbf{y}|\mathbf{b})}$
(8)

which is composed of the surprise plus a term capturing the distance from the
probability of external states given the blanket $p(\mathbf{y}|\mathbf{b})$ to
a variational model of the environment $q(\mathbf{y}|\boldsymbol{\theta})$
parametrized by $\boldsymbol{\theta}$.

This free energy constitutes a bound on the surprise, which is exact when the
variational distribution $q(\mathbf{y}|\boldsymbol{\theta})$ is equal to the
reference distribution $p(\mathbf{y}|\mathbf{b})$. The FEP often assumes that
$q$ is a normal distribution:

$\displaystyle q(\mathbf{y}|\boldsymbol{\theta})=$
$\displaystyle\mathcal{N}(\boldsymbol{\theta},\boldsymbol{\Sigma}^{\theta}),$
(9)

where $\boldsymbol{\theta}$ are the most likely states of the system and
$\boldsymbol{\Sigma}^{\theta}$ its covariance matrix.

A simple instance of Bayesian inference can be derived from the assumption
that the integral of the negative conditional log-probability (or surprise) in
the Free energy equation is approximately quadratic in the region near the
mode of the conditional density $\boldsymbol{\theta}$. This is called the
Laplace approximation [44] and allows approximating

$\displaystyle F(\boldsymbol{\theta},\mathbf{b})\approx$
$\displaystyle-\frac{n_{y}}{2}(1+\log(2\pi))-\frac{1}{2}\log|\boldsymbol{\Sigma}^{\theta}|$
$\displaystyle+\int
d\mathbf{y}q(\mathbf{y}|\boldsymbol{\theta})\left(\textfrak{J}(\boldsymbol{\theta},\mathbf{b})+\frac{1}{2}(\mathbf{y}-\boldsymbol{\theta})^{\intercal}\mathbf{H}_{yy}(\mathbf{y}-\boldsymbol{\theta})\right)$
$\displaystyle=$
$\displaystyle\textfrak{J}(\boldsymbol{\theta},\mathbf{b})+\frac{1}{2}\operatorname{Tr}\left[\boldsymbol{\Sigma}^{\theta}\mathbf{H}_{yy}\right]-\frac{n_{y}}{2}(1+\log(2\pi))-\frac{1}{2}\log|\Sigma^{\theta}|.$
(10)

where $\operatorname{Tr}$ is a trace operator,
$\mathbf{H}_{yy}=\nabla_{\mathbf{y}\mathbf{y}}\textfrak{J}(\boldsymbol{\theta},\mathbf{b})$
is the Hessian matrix respect to $\mathbf{y}$ of the marginal probability
distribution $p(\mathbf{y},\mathbf{b})$ at $\mathbf{y}=\boldsymbol{\theta}$.

In this scenario, finding the model $q(\mathbf{y}|\boldsymbol{\theta})$ that
minimizes the variational free energy is equivalent to approximating the
distribution $p(\mathbf{y}|\mathbf{b})$. From a gradient descent perspective,
if we are interested in adjusting parameters $\boldsymbol{\theta}$, this
minimization process results from following the negative gradient

$\displaystyle\nabla_{\boldsymbol{\theta}}F(\boldsymbol{\theta},\mathbf{b})=$
$\displaystyle\nabla_{\boldsymbol{\theta}}\textfrak{J}(\boldsymbol{\theta},\mathbf{b}).$
(11)

Under the FEP it is suggested that a system displaying a Markov blanket
minimizes the free energy functional by implementing a gradient descent scheme
referred to as recognition dynamics [7, 45, 46]. A literal interpretation of
this claim involves a dynamics of the variational parameter with the form

$\displaystyle\frac{\mathrm{d}\boldsymbol{\theta}_{t}}{\mathrm{d}t}=$
$\displaystyle-\boldsymbol{\gamma}\nabla_{\boldsymbol{\theta}}\textfrak{J}(\boldsymbol{\theta}_{t},\mathbf{b}_{t}),$
(12)

where $\boldsymbol{\gamma}$ is a matrix characterizing the rates of adjustment
of the variational parameters. This type of dynamics (or a discrete
counterpart) is usually proposed by active inference schemes (e.g. [46]).
However, in the most recent articles (e.g. [30]) this claim is not taken
literally, and the proponents of the FEP suggest that it is only the average
flows of the system (not the actual dynamics) that point in the direction of
the free energy gradient. Throughout this manuscript, we will explore both a
literal interpretation of a gradient descent on the free energy and the case
in which the free energy gradient is only connected with the average flows.

The FEP asserts that any system with a Markov blanket partition that reaches a
non-equilibrium steady-state (NESS) can be construed as performing an
elementary sort of Bayesian inference. This implies that the behaviour of a
system can be described by some variable that behaves as
$\boldsymbol{\theta}_{t}$ in Eq. 12, at least on average. Specifically, the
FEP proposes that the evolution of the statistics of internal states can be
described in terms of the variational free energy, given a blanket
$\mathbf{b}_{t}$.

Under Condition 3, given a blanket state at time $t$ ($\mathbf{b}_{t}$), the
statistics of internal and external states can be described independently. The
FEP proposes to describe change in internal and external states through
variables encoding the most likely222Choosing the $\mathrm{arg\,max}$ function
can imply problems in some cases (e.g. if it is non-differentiable), and in
some cases the expectation has been proposed as an alternative statistic.
However, in the case of Gaussian systems these two functions are equivalent,
so this does not affect the conclusions in this article. internal and external
states conditioned on the blanket

$\displaystyle\mathbf{m}_{x}(\mathbf{b}_{t})=$
$\displaystyle\mathrm{arg\,max}_{\mathbf{x}}p(\boldsymbol{x}_{t}|\mathbf{b}_{t}),$
(13) $\displaystyle\mathbf{m}_{y}(\mathbf{b}_{t})=$
$\displaystyle\mathrm{arg\,max}_{\mathbf{y}}p(\mathbf{y}_{t}|\mathbf{b}_{t}).$
(14)

Then, the average flow of the system (or marginal flow) conditioned on the
blanket $\mathbf{b}_{t}$ can be computed from the SDE decomposition in
Condition 1 as

$\displaystyle\left\langle
f_{y}(\mathbf{y}_{t},\mathbf{b}_{t})\right\rangle_{\mathbf{b}_{t}}=$
$\displaystyle\int
d\mathbf{x}_{t}d\mathbf{y}_{t}p(\mathbf{y}_{t},\mathbf{x}_{t}|\mathbf{b}_{t})f_{y}(\mathbf{y}_{t},\mathbf{b}_{t})$
$\displaystyle=$ $\displaystyle\int
d\mathbf{x}_{t}d\mathbf{y}_{t}p(\mathbf{y}_{t},\mathbf{x}_{t}|\mathbf{b}_{t})(\mathbf{Q}_{yz}-\boldsymbol{\Gamma}_{yz})\nabla_{\mathbf{z}}\textfrak{J}(\mathbf{z}_{t})$
$\displaystyle=$
$\displaystyle(\mathbf{Q}_{yy}-\boldsymbol{\Gamma}_{yy})\left\langle\nabla_{\mathbf{y}}\textfrak{J}(\mathbf{z}_{t})\right\rangle_{\mathbf{b}_{t}}$
$\displaystyle+\mathbf{Q}_{yb}\left\langle\nabla_{\mathbf{b}}\textfrak{J}(\mathbf{z}_{t})\right\rangle_{\mathbf{b}_{t}}+\mathbf{Q}_{yx}\left\langle\nabla_{\mathbf{x}}\textfrak{J}(\mathbf{z}_{t})\right\rangle_{\mathbf{b}_{t}}.$
(15)

The first term in this expression can be related to the gradient of the
surprise in Eq. 11 (see below). However, the second and third terms in the Eq.
15 preclude a straightforward connection between the average flow of the
system and the minimization of the variational free energy. The necessary step
in deriving the FEP is removing the solenoidal couplings between blocks of the
system, encoded in the matrix $\mathbf{Q}$, to remove these second and third
terms. Thus, in order to describe the equivalence between the dynamics of a
system and free energy minimization, the FEP assumes that

######  Assumption 1.

Solenoidal couplings between ‘blocks’ of states
($\mathbf{y},\mathbf{s},\mathbf{a},\mathbf{x}$) are precluded when a Markov
blanket emerges under sparse coupling [3].

$\displaystyle\mathbf{Q}=\begin{bmatrix}\mathbf{Q}_{yy}&&&\\\
&\mathbf{Q}_{ss}&&\\\ &&\mathbf{Q}_{aa}&\\\ &&&\mathbf{Q}_{xx}\end{bmatrix}.$
(16)

This leads to

$\displaystyle\left\langle
f_{y}(\mathbf{y}_{t},\mathbf{b}_{t})\right\rangle_{\mathbf{b}_{t}}=$
$\displaystyle(\mathbf{Q}_{yy}-\boldsymbol{\Gamma}_{yy})\left\langle\nabla_{\mathbf{y}}\textfrak{J}(\mathbf{z}_{t})\right\rangle_{\mathbf{b}_{t}}$
$\displaystyle\approx$
$\displaystyle(\mathbf{Q}_{yy}-\boldsymbol{\Gamma}_{yy})\nabla_{\mathbf{m}_{y}}\textfrak{J}(\mathbf{m}_{y}(\mathbf{b}_{t}),\mathbf{b}_{t}).$
(17)

Where the approximation is obtained neglecting couplings of order larger than
quadratic in the average of the surprise $\textfrak{J}(\mathbf{z}_{t})$, as
prescribed by the Laplace assumption, knowing that the flow of external states
is independent of internal states $\mathbf{x}$. The obtained expression is
proportional to the gradient in Eq. 11 for
$\boldsymbol{\theta}_{t}=\mathbf{m}_{y}(\mathbf{b}_{t})$, therefore pointing
to a direction minimizing the free energy, where the factor
$-(\mathbf{Q}_{yy}-\boldsymbol{\Gamma}_{yy})$ represents the rate of
adjustment of variational parameters ($\boldsymbol{\gamma}$ in Eq. 12).

In a recent work [30], it has been proposed that in the general case
$\mathbf{Q}$ is only block-diagonal for autonomous
($\\{\mathbf{a},\mathbf{x}\\}$) and ‘non-autonomous’ states
($\\{\mathbf{y},\mathbf{s}\\}$), allowing non-zero components
$\mathbf{Q}_{ys},\mathbf{Q}_{sy},\mathbf{Q}_{ax},\mathbf{Q}_{xa}$. However, it
is not clear (to our knowledge) how this case can be directly connected with
free energy minimization, and this should be perhaps clarified by future work.
In any case, the findings in the next sections of this article bring forward
similar problems considering one type of block-diagonal matrix or the other.

This step concludes the first move for deriving the FEP, by connecting the
average flow of a system with the gradient minimizing a variational free
energy functional.

###  2.2 Second move: Linking the average flow with the dynamics of the most
likely states

The second move for deriving the FEP involves connecting the average flow in
the system with its (averaged) dynamics. This second step starts by assuming a
mapping connecting the most likely internal and external states:

######  Assumption 2.

There is a smooth and differentiable function $\sigma$ that maps between the
most likely internal and external states given a blanket state,

$\displaystyle\mathbf{m}_{y}(\mathbf{b}_{t})=\sigma(\mathbf{m}_{x}(\mathbf{b}_{t})),$
(18)

and the gradient $\nabla_{\mathbf{m}_{x}}\sigma(\mathbf{m}_{x})$ is invertible
(i.e. $(\nabla_{\mathbf{m}_{x}}\sigma(\mathbf{m}_{x}))^{-1}$ exists)

A sufficient requisite for Assumption 2 is that the mapping from
$\mathbf{b}_{t}$ to $\mathbf{m}_{x}(\mathbf{b}_{t})$ is injective [4].

Once the mapping between internal and external states is defined, the next
step, as we anticipated, admits two possible interpretations. The first is an
interpretation in which the dynamics of the average states of the system
strictly follow a gradient descent on the free energy (in the form of Eq. 12,
e.g. [46]). A second interpretation relaxes this view to propose that free
energy minimization only takes place on average over counterfactual
trajectories (rather than directly). The distinction between the two (see
Assumptions 3∗ and 3∗∗) has been generally not discussed in detail, but it is
of great importance to evaluate the claims of the FEP.

The first interpretation proposes that the dynamics of the most likely states
can be described by the gradient on the free energy captured by the average
flow described by Eq. 17, which results in:

######  Assumption 3∗.

The evolution of the most likely external states is similar to their
conditional marginal flow given the blanket state

$\displaystyle\frac{\mathrm{d}\mathbf{m}_{y}(\mathbf{b}_{t})}{\mathrm{d}t}\approx$
$\displaystyle\left\langle
f_{y}(\mathbf{y}_{t},\mathbf{b}_{t})\right\rangle_{\mathbf{b}_{t}}.$ (19)

The star $*$ symbol in this assumption indicates that this assumption is in
general not explicitly stated, and that two competing interpretations are
possible. The first interpretation, in which Assumption 3∗ holds strictly, is
supported by verbal descriptions and some mathematical steps in [3] and [4]
(specifically the equivalents of Eq. 20 and 21 below333In these works, the
same symbol is used to represent
$\frac{\mathrm{d}\mathbf{m}_{y}(\mathbf{b}_{t})}{\mathrm{d}t}$ and
$\left\langle
f_{y}(\mathbf{y}_{t},\mathbf{b}_{t})\right\rangle_{\mathbf{b}_{t}}$ in Eq. 20
and 21 (Eq. 3.3 and 3.5 in [4], Eq. 8.22 and 8.23 in [3]), therefore making
the implicit assumption that they are approximately the same quantity. Also,
it should be noted that Eq. 22 can only be obtained under the combination of
Assumption 2 and Assumption 3∗, as it combines properties of the mapping
$\sigma$ (related to the most likely states) with the gradient of the free
energy (related to the average flows). and several verbal
descriptions444Particularly sentences about variable
$\dot{\boldsymbol{\eta}}(b)$ in these works, for which it is stated ‘the
dynamics of the internal mode as a gradient flow on the surprisal of the
external mode’ ([3], p. 97) ‘the rate of change of the most likely internal
states’ ([4]), p. 7.)). The alternative interpretation of this assumption
relaxes the equivalence between
$\frac{\mathrm{d}\mathbf{m}_{y}(\mathbf{b}_{t})}{\mathrm{d}t}$ and
$\left\langle
f_{y}(\mathbf{y}_{t},\mathbf{b}_{t})\right\rangle_{\mathbf{b}_{t}}$ (see
Assumption 3∗∗ below).

Taken together, these three assumptions let us derive a connection between the
evolution of the most likely states and the gradient of the variational free
energy:

$\displaystyle\frac{\mathrm{d}\mathbf{m}_{y}(\mathbf{b}_{t})}{\mathrm{d}t}=$
$\displaystyle(\mathbf{Q}_{yy}-\boldsymbol{\Gamma}_{yy})\nabla_{\mathbf{m}_{y}}\textfrak{J}(\mathbf{m}_{y}(\mathbf{b}_{t}),\mathbf{b}_{t}).$
(20)

This has an important implication, as it allows us to derive that
$\mathbf{m}_{y}(\mathbf{b}_{t})$ behaves as the variational inference
parameter $\boldsymbol{\theta}$ in Eq. 11. Therefore, the dynamics of a system
can be described as if performing variational inference.

This is possible because the conditional synchronisation manifold (Assumption
2) allows deriving a mapping between the evolution of the most likely states
through the chain rule:

$\displaystyle\mathbf{m}_{y}(\mathbf{b}_{t})=\sigma(\mathbf{m}_{x}(\mathbf{b}_{t}))\Rightarrow\frac{\mathrm{d}\mathbf{m}_{y}(\mathbf{b}_{t})}{\mathrm{d}t}=\nabla_{\mathbf{m}_{x}}\sigma\frac{\mathrm{d}\mathbf{m}_{x}(\mathbf{b}_{t})}{\mathrm{d}t}.$
(21)

Finally, under Assumption 3∗, the dynamics of the most likely internal states
can be described as minimizing the variational free energy about external
states

$\displaystyle\frac{\mathrm{d}\mathbf{m}_{x}(\mathbf{b}_{t})}{\mathrm{d}t}=$
$\displaystyle(\nabla_{\mathbf{m}_{x}}\sigma)^{-1}\frac{\mathrm{d}\mathbf{m}_{y}(\mathbf{b}_{t})}{\mathrm{d}t}$
$\displaystyle\approx$
$\displaystyle(\nabla_{\mathbf{m}_{x}}\sigma)^{-1}(\mathbf{Q}-\boldsymbol{\Gamma})(\nabla_{\mathbf{m}_{x}}\sigma)^{-1}\nabla_{\mathbf{m}_{x}}\textfrak{J}(\sigma(\mathbf{m}_{x}(\mathbf{b}_{t})),\mathbf{b}_{t}),$
(22)

which corresponds to following the gradient descent on the variational free
energy described by Eq. 12, now rewritten in terms of the most likely internal
states.

Here, we shall note that, in general, Assumption 3∗ does not hold for most
dynamical systems with stochastic fluctuations, as it equates the rate of
change of an average with the average of the rate of change. The proponents of
the FEP in recent works offer instead a more relaxed interpretation of the
gradient descent on the variational free energy. This interpretation proposes
that the FEP applies just to the marginal flows, and thus a system behaves ‘as
if’ performing Bayesian inference just on average. For example in [30] the
authors propose that ‘the interpretation in terms of Bayesian inference
emerges only in expectation –- or on average […] The classical example here is
the averaging of multiple responses to sensory perturbations, when
characterizing evoked responses in internal states’. This results in the
substitution of Assumption 3∗ by

######  Assumption 3∗∗.

If the conditional average flows follow the direction of a descending gradient
of a variational free energy, the behaviour of the states of the system can be
interpreted ‘as if’ they were, on average, performing a gradient descent or
minimizing a free energy functional.

In general, it is not easy to distinguish in the FEP literature when
Assumption 3∗ or Assumption 3∗∗ is considered, as verbal descriptions
sometimes refer to dynamics of most likely states and average flows
indistinctly. Note also that some important steps like the chain rule in Eq.
21 can only be derived for the interpretation promoted by Assumption 3∗.
Despite the interpretation, under these assumptions the proponents of the FEP
conclude that the dynamics of the most likely internal and external states can
be described as following a negative free energy gradient, which is equivalent
to stating that they evolve _as if_ performing a Bayesian inference. Moreover,
the FEP proposes that this can be extended to the dynamics of actions
$\mathbf{a}$ [4], deriving the principle of active inference [10, 12].

The two moves described here make important assumptions about the underlying
dynamical systems used to derive the theory. The first move assumes a very
specific statistical structure in which a gradient of the free energy
functional is directly connected to average flows in the systems, without
justifying to what extent it can be expected from the classes of systems
capturing properties from biological systems. The second move makes further
assumptions to justify that average flows (i.e. the average of the rate of
change) in the system are informative about its behaviour and dynamics (i.e.
the rate of change of the average), supporting an interpretation that
described the behaviour of a system _as if_ performing Bayesian inference. In
the next sections, we will see that the steps for deriving this interpretation
are not justified for the non-equilibrium linear systems studied in this
paper.

In the rest of the document, we will take Conditions 1 and 1 for granted and
explore the generality of the other conditions and assumptions in linear
stochastic systems with weak couplings. We will see that Assumption 2 hold
under special conditions and can be expected for a given class of systems. In
contrast, we show that Condition 3 and Assumption 1 only hold for very
specific sensorimotor loops and Assumption 3∗ and 3∗∗ do not hold in general,
threatening the viability of the FEP and its applicability to most living
systems.

##  3 Mathematical review of the FEP under linear stochastic dynamics

In order to explore the assumptions enumerated above in a class of tractable
non-equilibrium dynamical systems, we restrict our analysis to the class of
systems captured by a linear Langevin dynamics (which can be seen as an
approximation of an Ornstein-Uhlenbeck process [47]) defined by

$\frac{\mathrm{d}\mathbf{z}_{t}}{\mathrm{d}t}=\mathbf{J}(\mathbf{z}_{t}-\boldsymbol{\rho})+\boldsymbol{\omega}_{t},$
(23)

where $\mathbf{J}$ is an $n\times n$ invertible real matrix,
$\boldsymbol{\rho}$ is an $n$ dimensional real vector,
$\boldsymbol{\omega}_{t}$ is a standard n-dimensional Gaussian white noise
with a diagonal covariance matrix $2\boldsymbol{\Gamma}$. The linearity of the
process guarantees that the model will eventually result in a Gaussian
distribution.

The solution of this system (see B) in the non-equilibrium steady state (NESS)
takes the form of a multivariate Gaussian distribution
$\mathcal{N}(\boldsymbol{\rho},\boldsymbol{\Sigma}^{*})$ [47, 48] with
statistical moments:

$\displaystyle\lim_{t\to\infty}\mathbf{m}_{t}=$
$\displaystyle\boldsymbol{\rho},$ (24)
$\displaystyle\lim_{t\to\infty}\boldsymbol{\Sigma}_{t}=$
$\displaystyle\boldsymbol{\Sigma}^{*},\qquad\mathbf{J}\boldsymbol{\Sigma}^{*}+\boldsymbol{\Sigma}^{*}\mathbf{J}^{\intercal}+2\boldsymbol{\Gamma}=\mathbf{0},$
(25)

where $\boldsymbol{\Sigma}^{*}$ can be found numerically by solving the above
_continuous Lyapunov_ equation. If $\mathbf{J}$ is symmetric, the steady state
of the system is a state of equilibrium with
$\boldsymbol{\Sigma}^{*}=-\mathbf{J}^{-1}\boldsymbol{\Gamma}$. However, the
FEP focuses instead on NESS, which are more appropriate for describing living
systems.

For studying the NESS of the system, in this article we explore the case in
which non-diagonal couplings in $\mathbf{J}$ are small, although we explore
the effect of considering higher orders of this approximation. Thus, we define
the coupling matrix as, $\mathbf{J}=-\mathbf{I}+\mathbf{C}$, where
$\mathbf{C}^{2}$ are assumed to be small.

This leads to the following power series expansion,

$\displaystyle\boldsymbol{\Sigma}^{*}=\boldsymbol{\Gamma}+\frac{1}{2}(\mathbf{C}\boldsymbol{\Gamma}+\boldsymbol{\Gamma}\mathbf{C}^{\intercal})+\frac{1}{4}\left(\mathbf{C}^{2}\boldsymbol{\Gamma}+2\mathbf{C}\boldsymbol{\Gamma}\mathbf{C}^{\intercal}+\boldsymbol{\Gamma}(\mathbf{C}^{\intercal})^{2}\right)+\mathcal{O}(\mathbf{C}^{3}).$
(26)

To further simplify things, we also define a homogeneously distributed noise
$\boldsymbol{\Gamma}=\varsigma^{2}\mathbf{I}$, being $\varsigma$ a scalar
constant. The details of the derivation of the equations above are described
in B.

Figure 4: Sensorimotor loop structures. The figure shows the ‘canonical’ or
general sensorimotor structure proposed by the FEP [30] (A), a circular loop
where all connections between elements are asymmetric and unidirectional (B)
and a restricted sensorimotor loop where all interactions between blocks of
the system are symmetric, only allowing asymmetric connections within blocks
(B).

###  3.1 Can we expect the requirements for deriving the FEP in living
systems?

To explore the generality of the conditions required to connect the average
flow of a system to the gradient of the free energy, here we explore
Conditions 3 and Assumption 1 required to derive the first move of the FEP.

#### Condition 3: Markov blanket

First, the FEP requires the existence of a Markov blanket imposing a
conditional independence between internal and external states given a blanket
state (Condition 3, [3, 4]). We will see that not all linear systems meet this
condition, although it can be considered as an approximation in the case of
very weak-couplings.

For systems represented by Gaussian distributions, the Markov blanket
condition (Condition 3) is only met when the inverse of the covariance, the
Hessian matrix $\mathbf{H}=\left(\boldsymbol{\Sigma}^{*}\right)^{-1}$,
satisfies,

$\displaystyle\mathbf{H}_{yx}=\mathbf{H}_{xy}=\mathbf{0}.$ (27)

Thus this begs the questions: How common are Markov blankets? And when can we
expect to find them in living systems? The FEP generally proposes that the
theory holds when the Markov blanket condition is satisfied under a ‘canonical
flow constraint’ [30] that is defined in our linear system by
$\nabla_{\mathbf{z}}f(\mathbf{z})=\mathbf{J}$ structured as in Fig. 4.A. This
means that, mechanistically, external and sensory states do not depend on
internal states, and that action and internal states do not depend on external
states. This is representative of the kind of asymmetries we expect from a
sensorimotor loop in living systems. However, even in linear systems, it is
not easy to know under what conditions a system can satisfy both the canonical
flow constraint and possess a Markov blanket. As [34] points out, neither one
of these conditions is sufficient to guarantee the other.

In the case of weak couplings where $\mathbf{J}=-\mathbf{I}+\mathbf{C}$, and
homogenous noise $\boldsymbol{\Gamma}=\varsigma^{2}\mathbf{I}$, the covariance
of the system can be expanded as Eq. 26, and the inverse covariance (Hessian)
can be computed as a Neumann series (B),

$\displaystyle\mathbf{H}=$
$\displaystyle\left(\boldsymbol{\Sigma}^{*}\right)^{-1}=\varsigma^{-2}\mathbf{I}-\frac{\varsigma^{-2}}{2}(\mathbf{C}+\mathbf{C}^{\intercal})+\frac{\varsigma^{-2}}{4}\left(\mathbf{C}^{\intercal}\mathbf{C}-\mathbf{C}\mathbf{C}^{\intercal}\right)+\mathcal{O}(\mathbf{C}^{3}).$
(28)

Under the couplings determined by Condition 1 we can see that for a first
order approximation $\mathbf{H}_{yx}=\mathbf{0}+\mathcal{O}(\mathbf{C}^{2})$,
satisfying the Markov blanket condition.

For a second order weak coupling approximation we have,

$\displaystyle\mathbf{H}_{yx}=-\frac{\varsigma^{-2}}{4}\left(\mathbf{C}_{ys}\mathbf{C}_{xs}^{\intercal}+\mathbf{C}_{ya}\mathbf{C}_{xa}^{\intercal}\right)+\mathcal{O}(\mathbf{C}^{3}).$
(29)

Under the canonical flow constraints, relatively few systems will display an
exact Markov blanket, except for combinations of parameters that happen to
cancel the terms in the equation above. One exception is systems with weak
couplings under circular loops (Fig. 4.B) or systems with two layers of
blanket states (e.g. the system in Fig. 4.C). Note that this is because cycles
generating conditional couplings between $\mathbf{x},\mathbf{y}$ are of order
higher than 2. In general, these cases will not display a Markov blanket for
stronger couplings (see Eq. 77,78, except for perfectly symmetric couplings).
Thus, we can conclude that Markov blankets will emerge only for particular
combinations of parameters, as cycles in the system will in general introduce
couplings preventing their existence.

#### Assumption 1: solenoidal coupling

The FEP requires that the averaged marginal flow of external states
$\mathbf{y}$, given a blanket state $\mathbf{b}$, depends only on the
gradients of its marginal density. For this, it is a requirement that there is
no solenoidal coupling between external and other states (Assumption 1). We
will see in this section that most linear systems will not meet this
condition.

We can rewrite Eq. 25 to express the matrix of solenoidal couplings
$\mathbf{Q}$ as (see B.2 or [34])

$\displaystyle\mathbf{J}\mathbf{Q}+\mathbf{Q}\mathbf{J}^{\intercal}=\mathbf{J}\boldsymbol{\Gamma}-\boldsymbol{\Gamma}\mathbf{J}^{\intercal}.$
(30)

Again, the values of $\mathbf{Q}$ can be obtained by solving the corresponding
continuous Lyapunov equation. Assuming $\mathbf{J}=\mathbf{C}-\mathbf{I}$,
this matrix can be expressed as the power series

$\displaystyle\mathbf{Q}=\frac{1}{2}(\mathbf{C}\boldsymbol{\Gamma}-\boldsymbol{\Gamma}\mathbf{C}^{\intercal})+\frac{1}{4}\left(\mathbf{C}^{2}\boldsymbol{\Gamma}-\boldsymbol{\Gamma}(\mathbf{C}^{2})^{\intercal}\right)+\mathcal{O}(\mathbf{C}^{3}).$
(31)

The first order approximation under weak coupling
$\mathbf{J}=-\mathbf{I}+\mathbf{C}$ and
$\boldsymbol{\Gamma}=\varsigma^{2}\mathbf{I}$ and the canonical flow
constraints (Condition 1, Fig. 4.A) results in the solenoidal coupling matrix

$\displaystyle\mathbf{Q}=\frac{\varsigma^{2}}{2}\begin{bmatrix}\mathbf{C}_{yy}-\mathbf{C}_{yy}^{\intercal}&\mathbf{C}_{ys}-\mathbf{C}_{sy}^{\intercal}&\mathbf{C}_{ya}&\\\
\mathbf{C}_{sy}-\mathbf{C}_{ys}^{\intercal}&\mathbf{C}_{ss}-\mathbf{C}_{ss}^{\intercal}&\mathbf{C}_{sa}-\mathbf{C}_{as}^{\intercal}&-\mathbf{C}_{xs}^{\intercal}\\\
-\mathbf{C}_{ya}^{\intercal}&\mathbf{C}_{as}-\mathbf{C}_{sa}^{\intercal}&\mathbf{C}_{aa}-\mathbf{C}_{aa}^{\intercal}&\mathbf{C}_{ax}-\mathbf{C}_{xa}^{\intercal}\\\
&\mathbf{C}_{xs}&\mathbf{C}_{xa}-\mathbf{C}_{ax}^{\intercal}&\mathbf{C}_{xx}-\mathbf{C}_{xx}^{\intercal}\end{bmatrix}+\mathcal{O}(\mathbf{C}^{2}).$
(32)

Where Assumption 1 is not met for most parameter combinations.

In the best-case scenario, we can make many terms in the matrix above
disappear by making blocks symmetric when possible. This however does not
completely remove non-diagonal blocks, and leaves

$\displaystyle\mathbf{Q}=\frac{\varsigma^{2}}{2}\begin{bmatrix}\mathbf{C}_{yy}-\mathbf{C}_{yy}^{\intercal}&&\mathbf{C}_{ya}&\\\
&\mathbf{C}_{ss}-\mathbf{C}_{ss}^{\intercal}&&-\mathbf{C}_{xs}^{\intercal}\\\
-\mathbf{C}_{ya}^{\intercal}&&\mathbf{C}_{aa}-\mathbf{C}_{aa}^{\intercal}&\\\
&\mathbf{C}_{xs}&&\mathbf{C}_{xx}-\mathbf{C}_{xx}^{\intercal}\end{bmatrix}+\mathcal{O}(\mathbf{C}^{2}).$
(33)

Thus we observe that, even in a very weakly coupled system, the only way of
setting $\mathbf{Q}_{ya}\mathbf{Q}_{sx}$ and their antisymmetric counterparts
to zero is to effectively decouple some parts of the sensorimotor loop
completely that results in a ‘symmetric’ interaction loop with the form
displayed in Fig. 4.C). In this type of system, detailed balance is only
broken by interactions inside blocks
$\mathbf{y},\mathbf{s},\mathbf{a},\mathbf{x}$, as all couplings between blocks
are symmetric. Therefore, the system is driven out-of-equilibrium only by
internal tendencies of these blocks, not by their interactions between them.
This precludes for example the existence of asymmetric agent-environment
interactions, which may be crucial for many living processes as a mechanism
for generating qualitative differences between the inside and the outside of a
system, as well as for regulating exchanges of matter and energy with the
environment [35, 37, 38, 39].

For a second order approximation, we find the solenoidal coupling terms
between external and autonomous states:

$\displaystyle\mathbf{Q}_{ya}=$
$\displaystyle\frac{\varsigma^{2}}{4}\left(2\mathbf{C}_{ya}+\mathbf{C}_{yy}\mathbf{C}_{ya}+\mathbf{C}_{ys}\mathbf{C}_{sa}+\mathbf{C}_{ya}\mathbf{C}_{aa}-\mathbf{C}_{sy}^{\intercal}\mathbf{C}_{as}^{\intercal}\right)+\mathcal{O}(\mathbf{C}^{3}),$
(34) $\displaystyle\mathbf{Q}_{yx}=$
$\displaystyle\frac{\varsigma^{2}}{4}\left(\mathbf{C}_{ya}\mathbf{C}_{ax}-\mathbf{C}_{sy}^{\intercal}\mathbf{C}_{xs}^{\intercal}\right)+\mathcal{O}(\mathbf{C}^{3}),$
(35) $\displaystyle\mathbf{Q}_{sa}=$
$\displaystyle\frac{\varsigma^{2}}{4}\Big{(}2(\mathbf{C}_{sa}-\mathbf{C}_{sa}^{\intercal})+\mathbf{C}_{sy}\mathbf{C}_{ya}+\mathbf{C}_{ss}\mathbf{C}_{sa}+\mathbf{C}_{sa}\mathbf{C}_{aa}$
$\displaystyle-\mathbf{C}_{ss}^{\intercal}\mathbf{C}_{as}^{\intercal}-\mathbf{C}_{as}^{\intercal}\mathbf{C}_{aa}^{\intercal}-\mathbf{C}_{xs}^{\intercal}\mathbf{C}_{ax}^{\intercal}\Big{)}+\mathcal{O}(\mathbf{C}^{3}),$
(36) $\displaystyle\mathbf{Q}_{sx}=$
$\displaystyle\frac{\varsigma^{2}}{4}\left(-2\mathbf{C}_{xs}^{\intercal}+\mathbf{C}_{sa}\mathbf{C}_{ax}-\mathbf{C}_{ss}^{\intercal}\mathbf{C}_{xs}^{\intercal}-\mathbf{C}_{as}^{\intercal}\mathbf{C}_{xa}^{\intercal}-\mathbf{C}_{xs}^{\intercal}\mathbf{C}_{xx}^{\intercal}\right)+\mathcal{O}(\mathbf{C}^{3}).$
(37)

We can see that this matrix will in general only be block diagonal in very
specific cases, exemplifying how the presence of higher-order terms
complicates the uncoupling of solenoidal terms.

These results show how, in the case of linear, weakly-coupled systems, the
assumptions about the statistical structure of a system required by the FEP is
restricted to a very narrow space of parameters (i.e. values of $\mathbf{C}$).
In particular, removing solenoidal couplings between blocks of the system
requires a highly symmetric coupling structure (Fig. 4.C). This prevents the
application of the theory to many common structures found in biological
systems.

###  3.2 Can the FEP explain and describe the behaviour of living systems?

Above, we have shown that the conditions for connecting the conditional
average flow with the gradient of the free energy functional only hold for a
very narrow class of systems. Now we explore, in the cases where the previous
requirements are met, the connection between the conditional average flow and
the behaviour of the most likely states of the system. That is, are the
results of the FEP able to describe, explain or predict the dynamics of the
systems that conform to its assumptions? Here we will see that, while
Assumption 2 follows from a broad class of linear systems, one of the major
findings of our study is that Assumptions 3∗ and 3∗∗ present important issues
that prevent the results of the FEP to be good descriptions of the behaviour
of stochastic dynamical systems

#### Assumption 2: Mapping between internal and external statistics

A further requirement of the FEP is that a mapping $\sigma$ exists between
internal and external states (Assumption 2). We will see that in linear
systems, at least, there are many systems that meet this condition.

Given a NESS well characterized by Gaussian distributions over states
$p(\mathbf{z}_{t})=\mathcal{N}(\boldsymbol{\rho},\boldsymbol{\Sigma^{*}})$,
the conditional distributions given a particular blanket state
$\mathbf{b}=\\{\mathbf{s},\mathbf{a}\\}$ results in the most likely states
(Eq. 13 and 14) being described as

$\displaystyle\mathbf{m}_{x}(\mathbf{b}_{t})=$
$\displaystyle\boldsymbol{\rho}_{x}+\boldsymbol{\Sigma}^{*}_{xb}\left(\boldsymbol{\Sigma}^{*}_{bb}\right)^{g}(\mathbf{b}_{t}-\boldsymbol{\rho}_{b}),$
(38) $\displaystyle\mathbf{m}_{y}(\mathbf{b}_{t})=$
$\displaystyle\boldsymbol{\rho}_{y}+\boldsymbol{\Sigma}^{*}_{yb}\left(\boldsymbol{\Sigma}^{*}_{bb}\right)^{g}(\mathbf{b}_{t}-\boldsymbol{\rho}_{b}),$
(39)

where $\left(\boldsymbol{\Sigma}^{*}_{bb}\right)^{g}$ is a generalized inverse
matrix (which is equivalent to the inverse for nonsingular matrices). If
$\boldsymbol{\Sigma}^{*}_{xb}$ and $\boldsymbol{\Sigma}^{*}_{bb}$ are
nonsingular (this implies that $n_{x}\geq n_{b}$), we can derive the linear
mapping,

$\displaystyle\mathbf{m}_{y}(\mathbf{b}_{t})=$
$\displaystyle\boldsymbol{\rho}_{y}+\boldsymbol{\Sigma}^{*}_{yb}\left(\boldsymbol{\Sigma}^{*}_{xb}\right)^{-1}\left(\mathbf{m}_{x}(\mathbf{b}_{t})-\boldsymbol{\rho}_{x}\right)\equiv\sigma(\mathbf{m}_{x}(\mathbf{b}_{t})),$
(40)

which is invertible if $\boldsymbol{\Sigma}^{*}_{yb}$ is non-singular (for
this $n_{y}\geq n_{b}$ is required)

$\displaystyle\mathbf{m}_{x}(\mathbf{b}_{t})=$
$\displaystyle\mathbf{m}_{x}+\boldsymbol{\Sigma}^{*}_{xb}\left(\boldsymbol{\Sigma}^{*}_{yb}\right)^{-1}\left(\mathbf{m}_{y}(\mathbf{b}_{t})-\boldsymbol{\rho}_{y}\right)\equiv\sigma^{-1}(\mathbf{m}_{y}(\mathbf{b}_{t})).$
(41)

Yielding the evolution of the most likely states as,

$\displaystyle\frac{\mathrm{d}\mathbf{m}_{x}(\mathbf{b}_{t})}{\mathrm{d}t}=$
$\displaystyle\nabla_{\mathbf{m}_{y}}\sigma^{-1}(\mathbf{m}_{y}(\mathbf{b}_{t}))\frac{\mathrm{d}\mathbf{m}_{y}(\mathbf{b}_{t})}{\mathrm{d}t}=\boldsymbol{\Sigma}^{*}_{xb}\left(\boldsymbol{\Sigma}^{*}_{yb}\right)^{-1}\frac{\mathrm{d}\mathbf{m}_{y}(\mathbf{b}_{t})}{\mathrm{d}t}.$
(42)

This shows that, in general, a mapping between the most likely internal and
external states exists if the corresponding covariance submatrices are
invertible.

The FEP often states that the existence of this mapping is a consequence of a
Markov blanket [4], but we note that the existence of such a mapping is
independent of Conditions 1 and 3. The existence of a Markov blanket implies a
conditional covariance $\boldsymbol{\Sigma}_{yx}(\mathbf{b}_{t})=\mathbf{0}$,
but this does not affect the relation between internal and external states.
Thus, in linear systems any variable can potentially mediate an invertible
mapping, independently of a Markov blanket, if the corresponding covariance
submatrices are nonsingular555At the moment of writing this manuscript, we
have found that a similar result has been found independently in unpublished
work [49]..

#### Assumption 3∗: The dynamics and Bayesian inference

The final assumption for deriving the FEP is that the evolution of the most
likely external states is linked with the evolution of the marginal flow
conditioned on blanket states $\mathbf{b}_{t}$. In this section, we see that
this assumption does not hold in general.

The time derivative of $\mathbf{m}_{y}(\mathbf{b}_{t})$ yields,

$\displaystyle\frac{\mathrm{d}\mathbf{m}_{y}(\mathbf{b}_{t})}{\mathrm{d}t}=$
$\displaystyle\boldsymbol{\Sigma}^{*}_{yb}\left(\boldsymbol{\Sigma}^{*}_{bb}\right)^{g}\frac{\mathrm{d}\mathbf{b}_{t}}{\mathrm{d}t}=\boldsymbol{\Sigma}^{*}_{yb}\left(\boldsymbol{\Sigma}^{*}_{bb}\right)^{g}\left(\mathbf{J}_{bz}(\mathbf{z}_{t}-\boldsymbol{\rho})+\boldsymbol{\omega}_{b,t}\right).$
(43)

Conversely, the equation proposed by the FEP uses the marginal flow
(Assumption 3∗, see [4] or Appendix B in [3]). We represent the dynamics of a
variable driven by this marginal flow as

$\displaystyle\frac{\mathrm{d}\mathbf{\widetilde{m}}_{y}(\mathbf{b}_{t})}{\mathrm{d}t}\equiv$
$\displaystyle\left\langle
f_{y}(\mathbf{y}_{t},\mathbf{b}_{t})\right\rangle_{\mathbf{b}_{t}}=(\mathbf{Q}_{yz}-\boldsymbol{\Gamma}_{yz})\left\langle\nabla_{\mathbf{z}}\textfrak{J}(\mathbf{z}_{t})\right\rangle_{\mathbf{b}_{t}}$
$\displaystyle=$
$\displaystyle\mathbf{J}_{yy}(\mathbf{m}_{y}(\mathbf{b}_{t})-\boldsymbol{\rho}_{y})+\mathbf{J}_{yb}(\mathbf{b}_{t}-\boldsymbol{\rho}_{b}),$
(44)

where $\mathbf{\widetilde{m}}_{y}(\mathbf{b}_{t})$ captures the evolution of
the most likely dynamics in a system behaving in a way that strictly minimizes
the free energy. The difference between the behaviour of
$\mathbf{m}_{y}(\mathbf{b}_{t})$ and
$\mathbf{\widetilde{m}}_{y}(\mathbf{b}_{t})$ allows us to evaluate how
informative is the FEP about the behaviour of a system. If a system
approximately follows a free energy gradient, the behaviour of these two
variables will present some similarities in their evolution. This
approximation is crucial for the FEP as it is this marginal flow that is
connected with the gradient of the free energy (therefore equating dynamics
and inference, see Assumption 3∗). By using Eq. 39 this expression can be
rewritten simply as,

$\displaystyle\frac{\mathrm{d}\mathbf{\widetilde{m}}_{y}(\mathbf{b}_{t})}{\mathrm{d}t}=$
$\displaystyle\left(\mathbf{J}_{yy}\boldsymbol{\Sigma}^{*}_{yb}\left(\boldsymbol{\Sigma}^{*}_{bb}\right)^{g}+\mathbf{J}_{yb}\right)(\mathbf{b}_{t}-\boldsymbol{\rho}_{b}).$
(45)

We can thus see that these two equations represent quite different quantities
and that, in general, Assumption 3∗ will not hold for linear systems.

A | B
---|---
|
Figure 5: Divergence between system evolution and the marginal flow.
Comparison between the true evolution of the system
$\mathbf{m}_{y}(\mathbf{b}_{t})$ versus the marginal flow
$\mathbf{\widetilde{m}}_{y}(\mathbf{b}_{t})$ for (A) $\varsigma=0.1$ and
random couplings $C_{ij}=\pm 0.1$ with the connection scheme from Fig. 4.C and
(B) $\varsigma=10^{-3}$ and random couplings $C_{ij}=\pm 0.1$ with the
connection scheme from Fig. 4.B but with $\mathbf{C}_{yb}=\mathbf{0}$. In both
cases $n_{y}=n_{x}=2,n_{a}=n_{s}=1$. Note that different parameters (e.g.
reducing the noise) yield qualitatively similar results.

Furthermore, we can show that this equivalence is also incorrect even in the
case of weak couplings. Specifically, for the first order weak coupling
approximation with $\mathbf{J}=-\mathbf{I}+\mathbf{C}$ and
$\boldsymbol{\Gamma}=\varsigma^{2}\mathbf{I}$. Assuming
$\boldsymbol{\Sigma}^{*}_{bb}$ is nonsingular (i.e.
$\left(\boldsymbol{\Sigma}^{*}_{bb}\right)^{g}=\left(\boldsymbol{\Sigma}^{*}_{bb}\right)^{-1}$),
the weak coupling expansion of its inverse according to the Neumann series
(other expansions exist in the case of generalized inverses [50]) is,

$\displaystyle\left(\boldsymbol{\Sigma}^{*}_{bb}\right)^{-1}=$
$\displaystyle\varsigma^{-2}(\mathbf{I}-\frac{1}{2}(\mathbf{C}_{bb}+\mathbf{C}_{bb}^{\intercal}))+\mathcal{O}(\mathbf{C}^{2}).$
(46)

Using the weak-coupling approximation of $\boldsymbol{\Sigma}^{*}$ (Eq. 26),
Eq. 43 results in

$\displaystyle\frac{\mathrm{d}\mathbf{m}_{y}(\mathbf{b}_{t})}{\mathrm{d}t}=$
$\displaystyle\frac{1}{2}(\mathbf{C}_{yb}+\mathbf{C}_{by}^{\intercal})(\mathbf{I}-\frac{1}{2}(\mathbf{C}_{bb}+\mathbf{C}_{bb}^{\intercal}))\big{(}(-\mathbf{I}_{bz}+\mathbf{C}_{bz})(\mathbf{z}_{t}-\boldsymbol{\rho})$
$\displaystyle+\boldsymbol{\omega}_{b,t}\big{)}+\mathcal{O}(\mathbf{C}^{2})$
$\displaystyle=$
$\displaystyle-\frac{1}{2}(\mathbf{C}_{yb}+\mathbf{C}_{by}^{\intercal})(\mathbf{b}_{t}-\boldsymbol{\rho}_{b}-\boldsymbol{\omega}_{b,t})+\mathcal{O}(\mathbf{C}^{2}).$
(47)

In contrast, the marginal flow in Eq. 3.2 results in

$\displaystyle\frac{\mathrm{d}\mathbf{\widetilde{m}_{y}}(\mathbf{b}_{t})}{\mathrm{d}t}=$
$\displaystyle\Big{(}(-\mathbf{I}_{yy}+\mathbf{C}_{yy})\frac{1}{2}(\mathbf{C}_{yb}+\mathbf{C}_{by}^{\intercal})(\mathbf{I}-\frac{1}{2}(\mathbf{C}_{bb}+\mathbf{C}_{bb}^{\intercal}))$
$\displaystyle+\mathbf{C}_{yb}\Big{)}(\mathbf{b}_{t}-\boldsymbol{\rho}_{b})+\mathcal{O}(\mathbf{C}^{2})$
$\displaystyle=$
$\displaystyle-\frac{1}{2}(-\mathbf{C}_{yb}+\mathbf{C}_{by}^{\intercal})(\mathbf{b}_{t}-\boldsymbol{\rho}_{b})+\mathcal{O}(\mathbf{C}^{2}),$
(48)

which not only ignores the random fluctuations term, but also reverses the
sign of the term $\mathbf{C}_{yb}$. Note that
$\mathbf{C}_{yb}=\mathbf{C}_{by}^{\intercal}$ when we force solenoidal
uncoupling (Assumption 1, see Section 3.1), making
$\frac{\mathrm{d}\mathbf{\widetilde{m}_{y}}(\mathbf{b}_{t})}{\mathrm{d}t}=\mathbf{0}$
in that case (while
$\frac{\mathrm{d}\mathbf{m}_{y}(\mathbf{b}_{t})}{\mathrm{d}t}$ is non-zero).

For this weak coupling approximation, the only case in which the approximation
is valid is one in which $\varsigma=0$ (a deterministic system) and
$\mathbf{C}_{yb}=\mathbf{0}$ (a system where the agent just observes the
environment without affecting it). Similar expressions could be derived for
higher order approximations.

As an example of the dissimilarity of these quantities, in Fig. 5.A we display
$\mathbf{m}_{y}(\mathbf{b}_{t})$ and
$\mathbf{\widetilde{m}}_{y}(\mathbf{b}_{t})$ for arbitrary parameters
structured as the sensorimotor loop described in Fig. 4.C with random
couplings $C_{ij}=\pm 0.1$ (note that some weights are set to zero and others
forced to be symmetric) and $n_{y}=n_{x}=2,n_{a}=n_{s}=1,\varsigma=0.1$. We
see that $\mathbf{\widetilde{m}}_{y}(\mathbf{b}_{t})$ cannot capture the
structure in $\mathbf{m}_{y}(\mathbf{b}_{t})$. The reason behind this is that
the derivative of $\mathbf{\widetilde{m}_{y}}(\mathbf{b}_{t})$ accumulates an
error from fluctuations in the system that are not captured, displaying a
random walk behaviour that is absent in the real variable
$\mathbf{m}_{y}(\mathbf{b}_{t})$ . Similarly, in the most favourable case,
setting a very small noise $\varsigma=10^{-3}$ and
$\mathbf{C}_{yb}=\mathbf{0}$ (Figure 5.B), the situation is similar, as even
very small noise terms are sufficient for driving the two terms apart due to
the integration of random fluctuations over time.

#### Assumption 3∗∗: A way out? Problems with interpreting behaviour as
Bayesian inference only ‘on average’

The results of the previous section suggest that the FEP could be, in
practice, inapplicable for describing or explaining the behaviour of living
systems. Some of the most recent works on the FEP [30] try to avoid the
problems described above and state that the free energy principle describes
just marginal flows, that is, it describes the behaviour of a system just on
average over different trajectories. This could appear to circumvent some
issues presented in the previous section, as it implies substituting
Assumption 3∗ by a more relaxed interpretation of a gradient descent on the
free energy described by Assumption 3∗∗. However, under a close inspection, we
encounter that the situation is not improved by this claim. Assumption 3∗∗
entails two problems: 1) the mapping described by Assumption 2 no longer
connects internal and external flows, and 2) the conditional average flow
following the gradient of the free energy does not guarantee a gradient
descent in the effective behaviour. While the first problem can be solved, the
second has deep implications as in most cases the average flow does not
describe the true behaviour of a system in the presence of stochastic
fluctuations, therefore threatening to render the FEP inapplicable to describe
most living systems. In this last section, we briefly explore the validity of
the theory when applied only to the average flows of a system (and not the
evolution of the most likely states).

A first, practical problem implied by this claim is that if the FEP only
applies to marginal flows $\left\langle
f_{y}(\mathbf{y}_{t},\mathbf{b}_{t})\right\rangle_{\mathbf{b}_{t}},\left\langle
f_{x}(\mathbf{b}_{t},\mathbf{x}_{t})\right\rangle_{\mathbf{b}_{t}}$, then a
new mapping between flows is required. As we described in Eq. 21, a mapping
between the dynamics of the most likely states is derived from the chain rule
as the gradient of the mapping $\sigma$. If instead one interprets the FEP as
connecting the conditional average flow of external and internal states a new
mapping is required:

$\displaystyle\left\langle
f_{y}(\mathbf{y}_{t},\mathbf{b}_{t})\right\rangle_{\mathbf{b}_{t}}=\phi(\left\langle
f_{x}(\mathbf{b}_{t},\mathbf{x}_{t})\right\rangle_{\mathbf{b}_{t}}).$ (49)

In general, this mapping can take complicated forms, and often a unique
mapping will not exist. In linear systems, however, we can simplify the
marginal flows

$\displaystyle\left\langle
f_{y}(\mathbf{y}_{t},\mathbf{b}_{t})\right\rangle_{\mathbf{b}_{t}}=$
$\displaystyle\mathbf{J}_{yy}(\mathbf{m}_{y}(\mathbf{b}_{t})-\boldsymbol{\rho}_{y})+\mathbf{J}_{yb}(\mathbf{b}_{t}-\boldsymbol{\rho}_{b})$
$\displaystyle=$
$\displaystyle\left(\mathbf{J}_{yy}+(\boldsymbol{\Sigma}_{yb}^{*}(\boldsymbol{\Sigma}_{bb}^{*})^{g})^{-1}\right)(\mathbf{m}_{y}(\mathbf{b}_{t})-\boldsymbol{\rho}_{y})$
$\displaystyle=$
$\displaystyle\left(\mathbf{J}_{yy}+(\boldsymbol{\Sigma}_{yb}^{*}(\boldsymbol{\Sigma}_{bb}^{*})^{g})^{-1}\right)\boldsymbol{\Sigma}^{*}_{yb}\left(\boldsymbol{\Sigma}^{*}_{xb}\right)^{-1}\left(\mathbf{m}_{x}(\mathbf{b}_{t})-\boldsymbol{\rho}_{x}\right),$
(50) $\displaystyle\left\langle
f_{x}(\mathbf{b}_{t},\mathbf{x}_{t})\right\rangle_{\mathbf{b}_{t}}=$
$\displaystyle\mathbf{J}_{xx}(\mathbf{m}_{x}(\mathbf{b}_{t})-\boldsymbol{\rho}_{x})+\mathbf{J}_{xb}(\mathbf{b}_{t}-\boldsymbol{\rho}_{b})$
$\displaystyle=$
$\displaystyle\left(\mathbf{J}_{xx}+(\boldsymbol{\Sigma}_{xb}^{*}(\boldsymbol{\Sigma}_{bb}^{*})^{g})^{-1}\right)(\mathbf{m}_{x}(\mathbf{b}_{t})-\boldsymbol{\rho}_{x}).$
(51)

Which, if $\boldsymbol{\Sigma}_{bb}^{*}$ is invertible, yields

$\displaystyle\left\langle
f_{y}(\mathbf{y}_{t},\mathbf{b}_{t})\right\rangle_{\mathbf{b}_{t}}=$
$\displaystyle\left(\mathbf{J}_{yy}+\boldsymbol{\Sigma}_{bb}^{*}(\boldsymbol{\Sigma}_{yb}^{*})^{-1}\right)\boldsymbol{\Sigma}^{*}_{yb}\left(\boldsymbol{\Sigma}^{*}_{xb}\right)^{-1}$
$\displaystyle\cdot\left(\mathbf{J}_{xx}+\boldsymbol{\Sigma}_{bb}^{*}(\boldsymbol{\Sigma}_{xb}^{*})^{-1}\right)^{-1}\left\langle
f_{x}(\mathbf{b}_{t},\mathbf{x}_{t})\right\rangle_{\mathbf{b}_{t}}$
$\displaystyle\equiv$ $\displaystyle\phi(\left\langle
f_{x}(\mathbf{b}_{t},\mathbf{x}_{t})\right\rangle_{\mathbf{b}_{t}}),$ (52)

that results in the mapping $\phi$ being very different from
$\nabla_{\mathbf{m}_{x}}\sigma$. Thus, in general,

$\displaystyle\left\langle
f_{y}(\mathbf{y}_{t},\mathbf{b}_{t})\right\rangle_{\mathbf{b}_{t}}\neq\nabla_{\mathbf{m}_{x}}\sigma\left\langle
f_{x}(\mathbf{b}_{t},\mathbf{x}_{t})\right\rangle_{\mathbf{b}_{t}},$ (53)

thus contradicting [3, 4]. Thus, an interpretation of the FEP over conditional
average flows should replace Assumption 2 with a new mapping. This result, in
combination with a block-diagonal Q (Assumption 1) allows rewriting
$\left\langle
f_{y}(\mathbf{y}_{t},\mathbf{b}_{t})\right\rangle_{\mathbf{b}_{t}}$ and
$\left\langle
f_{x}(\mathbf{x}_{t},\mathbf{b}_{t})\right\rangle_{\mathbf{b}_{t}}$ pointing
in the direction of free energy minimization, mediated by a mapping $\phi$.

However, there is an important conceptual problem that remains even if a new
mapping is derived. The FEP relies on finding a variable that behaves as a
gradient descent on the free energy functional as described in Eq. 12. If
Assumption 3∗ is not met, then there is no variable in the system that behaves
following a gradient descent on the free energy (i.e. a variable
$\boldsymbol{\theta}_{t}$ with behaviour determined by
$\frac{\mathrm{d}\boldsymbol{\theta}_{t}}{\mathrm{d}t}=-\boldsymbol{\gamma}\nabla_{\boldsymbol{\theta}}F(\boldsymbol{\theta}_{t},\mathbf{b}_{t})$).
In this case, this assumption could be relaxed to Assumption 3∗∗ but, as we
observed in the results and simulations above, even when the average flow of
the system is related with the free energy (i.e.
$\left\langle\frac{\mathrm{d}\boldsymbol{\theta}_{t}}{\mathrm{d}t}\right\rangle_{\mathbf{b}_{t}}=-\boldsymbol{\gamma}\nabla_{\boldsymbol{\theta}}F(\boldsymbol{\theta}_{t},\mathbf{b}_{t})$)
this is not a good description of the behaviour of a system.

In this case, the claim of the FEP, if all other assumptions are hold, can be
described by Eq. 17, which relaxes the requirement of a strict gradient
descent (described by Eq. 20) and requires instead that the conditional
average flow is directed in the direction of the gradient of the free energy.
The problem with this relaxed gradient descent interpretation is that it does
not guarantee that a system effectively performs a gradient descent, not even
_on average_. In C we illustrate this issue in a simple bivariate linear
stochastic model with variables $y,b$. We observe that this particular model
presents a global attractor located at $y=b=0$ with solenoidal flows in the
form of a spiral flow (Fig. 2.A). In contrast, the conditional average flow
suggests a monotonic gradient ascent on $m_{y}^{2}(b)$, dismissing solenoidal
flows in the system and transforming an attracting flow into a repelling one.
This is a simple example showing how, in general, the conditional average
flows do not describe the behaviour of the system. Moreover, conditional
average flows can be misleading and not even a good approximation about the
average behaviour of a system, indicating a gradient ascent/descent on some
quantity when the behaviour of the system performs the opposite action. This
is also exemplified in recent work simulating linear systems, where the free
energy gradient only captures attracting tendencies at highly surprising
states, not capturing solenoidal flows nor behaviour near the NESS global
attractor [49, 43].

To summarize, our results show that, in general, even if a relaxed version of
the assumptions of the FEP holds, a system cannot be interpreted _as if_
performing Bayesian inference over external states. The reason behind this is
that the average flows $\left\langle
f_{y}(\mathbf{y}_{t},\mathbf{b}_{t})\right\rangle_{\mathbf{b}_{t}}$ and
$\left\langle
f_{x}(\mathbf{b}_{t},\mathbf{x}_{t})\right\rangle_{\mathbf{b}_{t}}$ do not
describe the behaviour of the system, as it can be easily seen from the
results in Fig. 5. Intuitively, substituting the true flow –
$f_{y}(\mathbf{y}_{t},\mathbf{b}_{t})$ – by an average flow fixing the blanket
state – $\left\langle
f_{y}(\mathbf{y}_{t},\mathbf{b}_{t})\right\rangle_{\mathbf{b}_{t}}$ –decouples
the trajectory of $\mathbf{y}$ from its previous state, which in most
dynamical systems will result in an impoverished description, not capturing
its real, history-dependent, behaviour.

##  4 Conclusion

The latest formulation of the free energy principle [3, 4] states that, in any
dynamical system equipped with a Markov blanket, the flow of internal states
can be construed as a gradient ascent on Bayesian model evidence. This
assumption rests on two crucial moves. The first move connects the system’s
average flow with a gradient on a variational free energy. This connection
relies on the existence of a Markov blanket and the emergence of a particular
statistical structure precluding solenoidal couplings. The second is the
interpretation that this relation between free energy and an average flow
results in systems behaving _as if_ performing variational inference over the
states of its environment. In this review, we have summarized crucial steps
required for this claim (Conditions 1–3 and Assumptions 1–3∗) and have shown
that several of these conditions cannot be met in general by linear, weakly-
coupled stochastic systems.

The first step compels a discussion about the generality of the FEP. That is,
if the principle requires a particular statistical structure, how general is
this structure, and how broadly can we expect it to be present among the class
of systems capturing the properties of living and cognitive processes? We
discover that, in the class of linear systems explored, the answer to this
question is that the statistical structure required by the FEP only arises in
a very narrow class of systems, requiring stringent conditions such as fully
symmetric agent-environment interactions that we cannot, in general, expect
from living systems [35, 37, 38, 39]. The generality of the FEP has been
questioned in the past due to conceptual issues [51, 33] or the existence of
counterexamples challenging the idea that perception-action interfaces, Markov
blankets and solenoidal decoupling follow from each other [34]. However, to
our knowledge, our study is the first that shows that the assumptions of the
FEP do not hold for a vast class of systems, namely, linear, weakly coupled
systems, except for the limited case of fully symmetric agent-environment
interaction.

This is concerning for two reasons. First, the FEP is designed for Gaussian
(i.e. in most cases linear) stochastic systems [4, 34]. Thus, our results
would imply that, as currently defined, the FEP cannot be fully implemented in
a broad set of systems from the class it was designed for. Second, one could
hope that the introduction of strong couplings or non-linearity allows some
specific systems to meet the required conditions. Some recent approximations
of the behaviour of chaotic oscillators point in this direction [42]. However,
this claim should be regarded with some scepticism, as the introduction of
stronger or higher-order couplings would result in additional terms in the
expansions explored in this article (see B). For most parameter settings,
systems with stronger couplings will not result in the independence relations
required for Condition 3 and Assumption 1. This will lead, in general, to a
more significant divergence between the evolution of the system and its
average flow, making it more unlikely that Assumptions 3∗ or 3∗∗ will hold. We
leave it to future work to explore the accuracy of this claim and investigate
whether the assumptions of the FEP can be met beyond the class of systems
explored here, e.g. displaying strong or nonlinear couplings.

The second step concerns how informative the FEP is about the behaviour of an
agent. The FEP justifies that any system can be described as if performing
variational inference through the existence of a conditional synchronization
manifold relating the direction of the free energy gradient and the evolution
of the most likely states of a system. We observe that this manifold can exist
in a broad class of systems. Nevertheless, assuming a strict gradient descent
interpretation (Assumption 3∗), it is problematic to connect the evolution of
the most likely states to the average flow of a system (and, in consequence,
to assume a system will behave as if minimizing the variational free energy).
The problem lies in implicitly relating the rate of change of the average
(expected) state as being described by the average flow (the expectation of
the rate of change) conditioned on a blanket state. If, instead, we consider a
more relaxed interpretation of the free energy gradient descent (i.e. just
taking place on average, Assumption 3∗∗), we encounter that new problems
arise. First, a new mapping between the flows of internal and external states
is required. Second, we observe that the average flow cannot, in general,
describe the true behaviour of a system. In sum, the FEP as it stands does not
do justice to the influence of the system’s trajectory in determining its
future behaviour. The reason behind this is that the gradient of the free
energy defined by the principle is computed for the average of an ensemble of
trajectories. Thus, even when the free energy gradient can be connected with
an average flow (which, as we have shown, happens under very specific
conditions), this is mainly uninformative about the behaviour of a system
subject to stochastic interactions. This is especially relevant for emergent
discussions about the compatibility of the FEP with enactive and autopoietic
theories of cognition (e.g. [52, 53, 54]). Specifically, enactive principles
stress the history-dependence of living systems, and this supposes a
fundamental incompatibility with the assumptions of the FEP [55]. In
particular, enactive views of cognition conceive sense-making as a process
emerging from the history of interactions of a system, which is invisible for
a gradient of the free energy described as an average flow.

The motivation behind the FEP aims to connect ideas from variational inference
with the dynamics of complex, self-organizing systems. This claim is
exceptionally appealing, as it could potentially allow applying the machinery
from Bayesian and information theoretical approaches to describe many systems
that are intractable in practice. However, by inspecting the theory and its
assumptions in the context of a broad class of analytically tractable models,
we discover that many of the steps required to derive the theory do not
straightforwardly follow or present significant conceptual problems that need
to be resolved. This finding illustrates the difficulties in developing a
theory of life and cognition over interdependent sets of mathematical
assumptions, and how testing these assumptions and their relations against
tractable models can help overcome these difficulties.

## Acknowledgements

We are grateful to Karl Friston, Lancelot Da Costa, Thomas Parr, Iñigo
Arandia-Romero and Hideaki Shimazaki for their constructive feedback and
comments on this manuscript. M.A. is thankful to Martin Biehl for helpful
discussions about references [3, 4, 30]. We are also grateful to Manuel
Baltieri for his open peer-review of a preprint of this manuscript.

M.A. was funded by the European Union’s Horizon 2020 research and innovation
programme under a Marie Skłodowska-Curie Action (grant agreement 892715) and
supported in part by the Basque Government project IT 1228-19 and the Spanish
Ministry of Science and Innovation project PID2019-104576GB-I00. C.L.B. is
supported by BBRSC grant BB/P022197/1.

## References

  * [1] K. Friston, P. Ao, Free energy, value, and attractors, Computational and mathematical methods in medicine 2012 (2012).
  * [2] K. Friston, Life as we know it, Journal of the Royal Society Interface 10 (86) (2013) 20130475.
  * [3] K. Friston, A free energy principle for a particular physics, arXiv preprint arXiv:1906.10184 (2019).
  * [4] T. Parr, L. Da Costa, K. Friston, Markov blankets, information geometry and stochastic thermodynamics, Philosophical Transactions of the Royal Society A 378 (2164) (2020) 20190159.
  * [5] J. Hohwy, The self-evidencing brain, Nous 50 (2) (2016) 259–285.
  * [6] A. Clark, Surfing uncertainty: Prediction, action, and the embodied mind, Oxford University Press, 2015.
  * [7] K. Friston, The free-energy principle: a rough guide to the brain?, Trends in cognitive sciences 13 (7) (2009) 293–301.
  * [8] K. Friston, J. Daunizeau, J. Kilner, S. J. Kiebel, Action and behavior: a free-energy formulation, Biological cybernetics 102 (3) (2010) 227–260.
  * [9] K. Friston, The free-energy principle: a unified brain theory?, Nature reviews neuroscience 11 (2) (2010) 127–138.
  * [10] K. Friston, J. Mattout, N. Trujillo-Barreto, J. Ashburner, W. Penny, Variational free energy and the laplace approximation, Neuroimage 34 (1) (2007) 220–234.
  * [11] K. Friston, A theory of cortical responses, Philosophical transactions of the Royal Society B: Biological sciences 360 (1456) (2005) 815–836.
  * [12] C. L. Buckley, C. S. Kim, S. McGregor, A. K. Seth, The free energy principle for action and perception: A mathematical review, Journal of Mathematical Psychology 81 (2017) 55–79.
  * [13] B. Millidge, A. Seth, C. L. Buckley, Predictive coding: a theoretical and experimental review, arXiv preprint arXiv:2107.12979 (2021).
  * [14] K. Friston, F. Rigoli, D. Ognibene, C. Mathys, T. Fitzgerald, G. Pezzulo, Active inference and epistemic value, Cognitive neuroscience 6 (4) (2015) 187–214.
  * [15] K. Friston, J. Daunizeau, S. J. Kiebel, Reinforcement learning or active inference?, PloS one 4 (7) (2009) e6421.
  * [16] K. Friston, T. FitzGerald, F. Rigoli, P. Schwartenbeck, G. Pezzulo, Active inference: a process theory, Neural computation 29 (1) (2017) 1–49.
  * [17] M. Cullen, B. Davey, K. Friston, R. J. Moran, Active inference in openai gym: A paradigm for computational investigations into psychiatric illness, Biological psychiatry: cognitive neuroscience and neuroimaging 3 (9) (2018) 809–818.
  * [18] B. Millidge, A. Tschantz, C. L. Buckley, Predictive coding approximates backprop along arbitrary computation graphs, arXiv preprint arXiv:2006.04182 (2020).
  * [19] J. Hohwy, A. Roepstorff, K. Friston, Predictive coding explains binocular rivalry: An epistemological review, Cognition 108 (3) (2008) 687–701.
  * [20] R. Kanai, Y. Komura, S. Shipp, K. Friston, Cerebral hierarchies: predictive processing, precision and the pulvinar, Philosophical Transactions of the Royal Society B: Biological Sciences 370 (1668) (2015) 20140169.
  * [21] A. Tschantz, B. Millidge, A. K. Seth, C. L. Buckley, Reinforcement learning through active inference, arXiv preprint arXiv:2002.12636 (2020).
  * [22] L. Da Costa, T. Parr, N. Sajid, S. Veselic, V. Neacsu, K. Friston, Active inference on discrete state-spaces: a synthesis, Journal of Mathematical Psychology 99 (2020) 102447.
  * [23] B. Millidge, Deep active inference as variational policy gradients, Journal of Mathematical Psychology 96 (2020) 102348.
  * [24] T. H. FitzGerald, P. Schwartenbeck, M. Moutoussis, R. J. Dolan, K. Friston, Active inference, evidence accumulation, and the urn task, Neural computation 27 (2) (2015) 306–328.
  * [25] T. Parr, D. Markovic, S. J. Kiebel, K. Friston, Neuronal message passing using mean-field, bethe, and marginal approximations, Scientific reports 9 (1) (2019) 1–18.
  * [26] D. Kappel, C. Tetzlaff, A synapse-centric account of the free energy principle, arXiv preprint arXiv:2103.12649 (2021).
  * [27] A. Tschantz, A. K. Seth, C. L. Buckley, Learning action-oriented models through active inference, PLoS computational biology 16 (4) (2020) e1007805.
  * [28] P. Calvo, K. Friston, Predicting green: really radical (plant) predictive processing, Journal of the Royal Society Interface 14 (131) (2017) 20170096.
  * [29] K. Friston, J. Kilner, L. Harrison, A free energy principle for the brain, Journal of physiology-Paris 100 (1-3) (2006) 70–87.
  * [30] K. Friston, L. Da Costa, T. Parr, Some interesting observations on the free energy principle, Entropy 23 (8) (2021) 1076.
  * [31] J. Pearl, Probabilistic reasoning in intelligent systems: networks of plausible inference, Elsevier, 1988.
  * [32] T. S. Richardson, P. Spirtes, et al., Automated discovery of linear feedback models, Carnegie Mellon, 1996.
  * [33] J. Bruineberg, K. Dolega, J. Dewhurst, M. Baltieri, The emperor’s new markov blankets, Behavioral and Brain Sciences (2021) 1–63.
  * [34] M. Biehl, F. A. Pollock, R. Kanai, A technical critique of some parts of the free energy principle, Entropy 23 (3) (2021) 293.
  * [35] J. E. Rothman, J. Lenard, Membrane asymmetry, Science 195 (4280) (1977) 743–753.
  * [36] G. Buzsaki, Rhythms of the Brain, Oxford University Press, 2006.
  * [37] B. Fadeel, D. Xue, The ins and outs of phospholipid asymmetry in the plasma membrane: roles in health and disease, Critical reviews in biochemistry and molecular biology 44 (5) (2009) 264–277.
  * [38] X. E. Barandiaran, E. Di Paolo, M. Rohde, Defining agency: Individuality, normativity, asymmetry, and spatio-temporality in action, Adaptive Behavior 17 (5) (2009) 367–386.
  * [39] K. Ruiz-Mirazo, A. Moreno, Basic autonomy as a fundamental step in the synthesis of life, Artificial life 10 (3) (2004) 235–259.
  * [40] C. Kwon, P. Ao, Nonequilibrium steady state of a stochastic system driven by a nonlinear drift force, Physical Review E 84 (6) (2011) 061106.
  * [41] R. Yuan, Y. Tang, P. Ao, Sde decomposition and a-type stochastic interpretation in nonequilibrium processes, Frontiers of Physics 12 (6) (2017) 1–9.
  * [42] K. Friston, C. Heins, K. Ueltzhöffer, L. Da Costa, T. Parr, Stochastic chaos and markov blankets, Entropy 23 (9) (2021) 1220.
  * [43] T. Parr, L. Da Costa, C. Heins, M. J. D. Ramstead, K. Friston, Memory and markov blankets, Entropy 23 (9) (2021) 1105.
  * [44] N. Reid, Approximate likelihoods, Proceedings of the ICIAM (2015).
  * [45] C. S. Kim, Recognition dynamics in the brain under the free energy principle, Neural computation 30 (10) (2018) 2616–2659.
  * [46] L. Da Costa, T. Parr, B. Sengupta, K. Friston, Neural dynamics under active inference: Plausibility and efficiency of information processing, Entropy 23 (4) (2021) 454.
  * [47] P. Vatiwutipong, N. Phewchean, Alternative way to derive the distribution of the multivariate ornstein–uhlenbeck process, Advances in Difference Equations 2019 (1) (2019) 1–7.
  * [48] C. Godrèche, J.-M. Luck, Characterising the nonequilibrium stationary states of ornstein–uhlenbeck processes, Journal of Physics A: Mathematical and Theoretical 52 (3) (2018) 035002.
  * [49] L. Da Costa, K. Friston, C. Heins, G. A. Pavliotis, Bayesian mechanics for stationary processes, arXiv preprint arXiv:2106.13830 (2021).
  * [50] J. J. Climent, N. Thome, Y. Wei, A geometrical approach on generalized inverses by neumann-type series, Linear algebra and its applications 332 (2001) 533–540.
  * [51] V. Raja, D. Valluri, E. Baggs, A. Chemero, M. L. Anderson, The markov blanket trick: On the scope of the free energy principle and active inference, Physics of Life Reviews (2021).
  * [52] J. Bruineberg, J. Kiverstein, E. Rietveld, The anticipating brain is not a scientist: the free-energy principle from an ecological-enactive perspective, Synthese 195 (6) (2018) 2417–2444.
  * [53] A. Constant, A. Clark, K. J. Friston, Representation wars: Enacting an armistice through active inference, Frontiers in Psychology 11 (2021) 3798.
  * [54] M. J. Ramstead, M. D. Kirchhoff, A. Constant, K. J. Friston, Multiscale integration: beyond internalism and externalism, Synthese 198 (1) (2021) 41–70.
  * [55] E. Di Paolo, E. Thompson, R. D. Beer, Laying down a forking path: Incompatibilities between enaction and the free energy principle, PsyArXiv (2021).

##  Appendix A Mathematical definitions of main concepts of the free energy
principle

Perception-action partition. The FEP assumes that the system can be decomposed
into external, sensory, active and internal states,
$\mathbf{z}=\\{\mathbf{y},\mathbf{s},\mathbf{a},\mathbf{x}\\}$, configured as
a perception-action loop reflecting an interface mediating between
‘autonomous’ states (active and internal states $\\{\mathbf{a},\mathbf{x}\\}$)
and ‘non-autonomous’ states (external and sensory states
$\\{\mathbf{y},\mathbf{s}\\}$). This leads to describing the evolution of the
system as

$\displaystyle\frac{\mathrm{d}\mathbf{z}_{t}}{\mathrm{d}t}=$ $\displaystyle
f(\mathbf{z}_{t})+\boldsymbol{\omega}_{t},$ (54) $\displaystyle
f(\mathbf{z}_{t})=$
$\displaystyle\begin{bmatrix}f_{y}(\mathbf{y}_{t},\mathbf{s}_{t},\mathbf{a}_{t})\\\
f_{s}(\mathbf{y}_{t},\mathbf{s}_{t},\mathbf{a}_{t})\\\
f_{a}(\mathbf{s}_{t},\mathbf{a}_{t},\mathbf{x}_{t})\\\
f_{x}(\mathbf{s}_{t},\mathbf{a}_{t},\mathbf{x}_{t})\end{bmatrix}.$ (55)

NESS and solenoidal flows. The FEP assumes that the system will reach a non-
equilibrium steady state described by the probability density function
$p(\mathbf{z}_{t})$, which can be described using a SDE decomposition that
separates the flow into dissipative
($-\boldsymbol{\Gamma}\nabla_{\mathbf{z}}\textfrak{J}(\mathbf{z}_{t})$) and
solenoidal ($\mathbf{Q}\nabla_{\mathbf{z}}\textfrak{J}(\mathbf{z}_{t})$)
components

$\displaystyle f(\mathbf{z}_{t})=$
$\displaystyle(\mathbf{Q}-\boldsymbol{\Gamma})\nabla_{\mathbf{z}}\textfrak{J}(\mathbf{z}_{t}),$
(56) $\displaystyle\textfrak{J}(\mathbf{z}_{t})=$ $\displaystyle-\log
p(\mathbf{z}_{t}),$ (57)

The FEP makes some requirements about the solenoidal flow matrix $\mathbf{Q}$,
as it assumes it is state-independent (i.e. $\mathbf{Q}$ does not change with
$\mathbf{z}$) and that it is sparse in the sense that couplings between some
states (e.g. internal and external) are zero (see 1).

Markov blanket. The steady state distribution of the system is described in
terms of a Markov blanket, where internal/external states are independent when
conditioned on its blanket states.

$\displaystyle
p(\mathbf{y}_{t},\mathbf{x}_{t}|\mathbf{b}_{t})=p(\mathbf{y}_{t}|\mathbf{b}_{t})p(\boldsymbol{x}_{t}|\mathbf{b}_{t}).$
(58)

Conditional synchronisation manifold. The FEP assumes that there is a smooth
and differentiable function, $\sigma$, which maps between the most likely
internal and external states given a blanket state,

$\displaystyle\mathbf{m}_{y}(\mathbf{b}_{t})=\sigma(\mathbf{m}_{x}(\mathbf{b}_{t})),$
(59)

and the gradient $\nabla_{\mathbf{m}_{x}}\sigma(\mathbf{m}_{x})$ is invertible
(i.e. $(\nabla_{\mathbf{m}_{x}}\sigma(\mathbf{m}_{x}))^{-1}$ exists). The FEP
generally refers to this mapping as a _conditional synchronisation manifold_ ,
and proposes that its existence allows to characterise the relationship
between (maximum a posteriori) internal and external states in terms of
internal states ‘sensing’ or ‘tracking’ external states through the Markov
blanket [3].

Rate of change of the average. By virtue of the conditional synchronisation
manifold, the rate of change of the average internal and external states of
the system are connected by the gradients of the mapping function.

$\displaystyle\mathbf{m}_{y}(\mathbf{b}_{t})=\sigma(\mathbf{m}_{x}(\mathbf{b}_{t}))\Rightarrow\frac{\mathrm{d}\mathbf{m}_{y}(\mathbf{b}_{t})}{\mathrm{d}t}=\nabla_{\mathbf{m}_{x}}\sigma\frac{\mathrm{d}\mathbf{m}_{x}(\mathbf{b}_{t})}{\mathrm{d}t}.$
(60)

Variational free energy. A system performing Bayesian inference tries to
minimize the surprise of observed states ($-\log p(\mathbf{b})$), according to
an internal model. However, a system cannot access this surprise value without
complete knowledge of its environment, Bayesian inference prescribes to use a
lower bound of this surprise described by the variational free energy
$F(\boldsymbol{\theta},\mathbf{b})$:

$\displaystyle-\log p(\mathbf{b})<F(\boldsymbol{\theta},\mathbf{b})=$
$\displaystyle-\log
p(\mathbf{b})+D_{KL}(q(\mathbf{y}|\boldsymbol{\theta})||p(\mathbf{y}|\mathbf{b})),$
(61)

which is composed of the surprise plus a term capturing the distance from the
probability of external states given the blanket $p(\mathbf{y}|\mathbf{b})$ to
an internal model of the world $q(\mathbf{y}|\boldsymbol{\theta})$
parametrized by $\boldsymbol{\theta}$.

Average flow and the free energy gradient. The FEP proposes that the evolution
of internal and external states of a system can be described, under a Gaussian
approximation, as a gradient of the free energy of the system respect to its
sufficient statistics.

$\displaystyle\left\langle
f_{y}(\mathbf{y}_{t},\mathbf{b}_{t})\right\rangle_{\mathbf{b}_{t}}=$
$\displaystyle(\mathbf{Q}_{yy}-\boldsymbol{\Gamma}_{yy})\left\langle\nabla_{\mathbf{y}}\textfrak{J}(\mathbf{z}_{t})\right\rangle_{\mathbf{b}_{t}}$
$\displaystyle\approx$
$\displaystyle(\mathbf{Q}_{yy}-\boldsymbol{\Gamma}_{yy})\nabla_{\mathbf{m}_{y}}F(\boldsymbol{\theta},\mathbf{b}).$
(62)

Then, the conditional synchronization manifold is proposed to connect the
internal and external average flows (i.e. the gradients
$\nabla_{\mathbf{m}_{y}}$ and $\nabla_{\mathbf{m}_{x}}$), suggesting that
internal states are behaving as if trying to perform Bayesian inference over
external states.

##  Appendix B Solution of the linear Langevin dynamics.

We start with the Ornstein-Uhlenbeck process

$\mathrm{d}\mathbf{z}_{t}=\mathbf{J}(\mathbf{z}_{t}-\boldsymbol{\rho})\mathrm{d}t+\mathrm{d}\mathbf{w}_{t},$
(63)

which can be approximated by the equivalent Langevin dynamics

$\frac{\mathrm{d}\mathbf{z}_{t}}{\mathrm{d}t}=\mathbf{J}(\mathbf{z}_{t}-\boldsymbol{\rho})+\boldsymbol{\omega}_{t},$
(64)

where $\mathbf{J}$ is an $n\times n$ invertible real matrix,
$\boldsymbol{\rho}$ is an $n$ dimensional real vector, $\mathbf{w}_{t}$ is a
standard n-dimensional Wiener process, $\boldsymbol{\omega}_{t}$ is a standard
n-dimensional Gaussian white noise with covariance matrix
$2\boldsymbol{\Gamma}$.

The model can be solved using standard methods for systems of differential
equations

$\displaystyle\mathrm{e}^{-\mathbf{J}t}\frac{\mathrm{d}\mathbf{z}_{t}}{dt}-\mathrm{e}^{-\mathbf{J}t}\mathbf{J}(\mathbf{z}_{t}-\boldsymbol{\rho})=$
$\displaystyle\mathrm{e}^{-\mathbf{J}t}\boldsymbol{\omega}_{t},$ (65)
$\displaystyle\frac{\mathrm{d}}{dt}(\mathrm{e}^{-\mathbf{J}t}(\mathbf{z}_{t}-\boldsymbol{\rho}))=$
$\displaystyle\mathrm{e}^{-\mathbf{J}t}\boldsymbol{\omega}_{t},$ (66)
$\displaystyle\mathbf{z}_{t}=$
$\displaystyle\mathrm{e}^{\mathbf{J}t}(\mathbf{z}_{0}-\boldsymbol{\rho})+\boldsymbol{\rho}+\int_{0}^{t}\mathrm{d}t^{\prime}\mathrm{e}^{-\mathbf{J}(t^{\prime}-t)}\boldsymbol{\omega}_{t^{\prime}}.$
(67)

This solution of the system [47, 48] yields statistical moments

$\displaystyle\mathbf{m}_{t}=$ $\displaystyle\int
d\boldsymbol{\omega}p(\boldsymbol{\omega})\mathbf{z}_{t}$ $\displaystyle=$
$\displaystyle\mathrm{e}^{\mathbf{J}t}\mathbf{z}_{0}+(\mathbf{I}-\mathrm{e}^{\mathbf{J}t})\boldsymbol{\rho},$
(68) $\displaystyle\boldsymbol{\Sigma}_{t}=$ $\displaystyle\int
d\boldsymbol{\omega}p(\boldsymbol{\omega})(\mathbf{z}_{t}-\mathbf{m}_{t})(\mathbf{z}_{t}-\mathbf{m}_{t})^{\intercal}$
$\displaystyle=$ $\displaystyle\int
d\boldsymbol{\omega}p(\boldsymbol{\omega})\int_{0}^{t}\mathrm{d}t^{\prime}\int_{0}^{t}\mathrm{d}t^{\prime\prime}\mathrm{e}^{\mathbf{J}(t-t^{\prime})}\boldsymbol{\omega}_{t^{\prime}}\boldsymbol{\omega}^{\intercal}_{t^{\prime\prime}}\mathrm{e}^{\mathbf{J}^{\intercal}(t-t^{\prime\prime})}$
$\displaystyle=$ $\displaystyle
2\int_{0}^{t}\mathrm{d}t^{\prime}\mathrm{e}^{-\mathbf{J}(t^{\prime}-t)}\boldsymbol{\Gamma}\mathrm{e}^{-\mathbf{J}^{\intercal}(t^{\prime}-t)}.$
(69)

These equations are hard to solve analytically. However, we can find the
differential equations that result in equivalent integrals, obtaining the time
evolution of statistical moments

$\displaystyle\frac{\mathrm{d}\mathbf{m}_{t}}{dt}=$
$\displaystyle\mathrm{e}^{\mathbf{J}t}\mathbf{J}(\mathbf{z}_{0}-\boldsymbol{\rho})=\mathbf{J}(\mathbf{m}_{t}-\boldsymbol{\rho}),$
(70) $\displaystyle\frac{\mathrm{d}\boldsymbol{\Sigma}_{zz,t}}{dt}=$
$\displaystyle
2\boldsymbol{\Gamma}+\boldsymbol{\Sigma}_{zz,t}\mathbf{J}^{\intercal}+\mathbf{J}\boldsymbol{\Sigma}_{zz,t}.$
(71)

At equilibrium or at the NESS, when the solution is unique, the system
stabilizes to the values that make the derivatives equal to zero

$\displaystyle\lim_{t\to\infty}\mathbf{m}_{t}=$
$\displaystyle\boldsymbol{\rho},$ (72)
$\displaystyle\lim_{t\to\infty}\boldsymbol{\Sigma}_{zz,t}=$
$\displaystyle\boldsymbol{\Sigma}^{*},\qquad\mathbf{J}\boldsymbol{\Sigma}^{*}+\boldsymbol{\Sigma}^{*}\mathbf{J}^{\intercal}+2\boldsymbol{\Gamma}=0.$
(73)

where $\boldsymbol{\Sigma}^{*}$ can be found numerically solving the above
_continuous Lyapunov_ equation (a particular case of a _Sylvester_ equation).
If $\mathbf{J}$ is symmetric, the steady state of the system is a state of
equilibrium with
$\boldsymbol{\Sigma}^{*}=-\mathbf{J}^{-1}\boldsymbol{\Gamma}$. However, the
FEP focuses instead on NESS, which are more appropriate for describing living
systems.

###  B.1 Weak coupling approximation

In order to study the solution at the NESS, we assume weak couplings of the
form $\mathbf{J}=-\mathbf{I}+\mathbf{C}$, with $\mathbf{C}^{2}$ being small.
We derive

$\displaystyle 2\boldsymbol{\Gamma}=$
$\displaystyle-\mathbf{J}\boldsymbol{\Sigma}^{*}-\boldsymbol{\Sigma}^{*}\mathbf{J}^{\intercal}=2\boldsymbol{\Sigma}^{*}-\mathbf{C}\boldsymbol{\Sigma}^{*}-\boldsymbol{\Sigma}^{*}\mathbf{C}^{\intercal},$
(74) $\displaystyle\boldsymbol{\Sigma}^{*}=$
$\displaystyle\boldsymbol{\Gamma}+\frac{1}{2}(\mathbf{C}\boldsymbol{\Sigma}^{*}+\boldsymbol{\Sigma}^{*}\mathbf{C}^{\intercal}).$
(75)

Recursively substituting $\boldsymbol{\Sigma}^{*}$ in the right-hand side with
$\boldsymbol{\Gamma}+\frac{1}{2}\mathbf{C}\boldsymbol{\Sigma}^{*}+\frac{1}{2}\boldsymbol{\Sigma}^{*}\mathbf{C}^{\intercal}$,
we obtain the time series expansion

$\displaystyle\boldsymbol{\Sigma}^{*}=$
$\displaystyle\boldsymbol{\Gamma}+\frac{1}{2}(\mathbf{C}\boldsymbol{\Gamma}+\boldsymbol{\Gamma}\mathbf{C}^{\intercal})+\frac{1}{4}\left(\mathbf{C}^{2}\boldsymbol{\Gamma}+2\mathbf{C}\boldsymbol{\Gamma}\mathbf{C}^{\intercal}+\boldsymbol{\Gamma}(\mathbf{C}^{\intercal})^{2}\right)$
$\displaystyle+\frac{1}{8}\left(\mathbf{C}^{3}\boldsymbol{\Gamma}+3\mathbf{C}^{2}\boldsymbol{\Gamma}\mathbf{C}^{\intercal}+3\mathbf{C}\boldsymbol{\Gamma}(\mathbf{C}^{\intercal})^{2}+\boldsymbol{\Gamma}(\mathbf{C}^{\intercal})^{3}\right)+\mathcal{O}(\mathbf{C}^{4}).$
(76)

In the equilibrium case of symmetric couplings, the Hessian is trivially
$\mathbf{H}=\mathbf{J}\boldsymbol{\Gamma}^{-1}$. For systems at a NESS, for
$\mathbf{J}=-\mathbf{I}+\mathbf{C}$ and under a uniform noise
$\boldsymbol{\Gamma}=\varsigma^{-2}\mathbf{I}$, the inverse covariance
(Hessian) can be computed as a Neumann series,

$\displaystyle\mathbf{H}=$
$\displaystyle\left(\boldsymbol{\Sigma}^{*}\right)^{-1}=\varsigma^{-2}\sum_{k=0}^{\infty}\Big{(}-\frac{1}{2}(\mathbf{C}+\mathbf{C}^{\intercal})-\frac{1}{4}\left(\mathbf{C}^{2}+2\mathbf{C}\mathbf{C}^{\intercal}+(\mathbf{C}^{\intercal})^{2}\right)$
$\displaystyle-\frac{1}{8}\left(\mathbf{C}^{3}+3\mathbf{C}^{2}\mathbf{C}^{\intercal}+3\mathbf{C}(\mathbf{C}^{\intercal})^{2}+(\mathbf{C}^{\intercal})^{3}\right)+\mathcal{O}(\mathbf{C}^{4})\Big{)}^{k}$
$\displaystyle=$
$\displaystyle\varsigma^{-2}\mathbf{I}-\frac{\varsigma^{-2}}{2}(\mathbf{C}+\mathbf{C}^{\intercal})-\frac{\varsigma^{-2}}{4}\left(\mathbf{C}^{2}+2\mathbf{C}\mathbf{C}^{\intercal}+(\mathbf{C}^{\intercal})^{2}\right)$
$\displaystyle+\frac{\varsigma^{-2}}{4}\left(\mathbf{C}^{2}+\mathbf{C}\mathbf{C}^{\intercal}+\mathbf{C}^{\intercal}\mathbf{C}+(\mathbf{C}^{\intercal})^{2}\right)$
$\displaystyle-\frac{\varsigma^{-2}}{8}\left(\mathbf{C}^{3}+3\mathbf{C}^{2}\mathbf{C}^{\intercal}+3\mathbf{C}(\mathbf{C}^{\intercal})^{2}+(\mathbf{C}^{\intercal})^{3}\right)$
$\displaystyle-\frac{\varsigma^{-2}}{8}\big{(}\mathbf{C}^{3}+\mathbf{C}^{2}\mathbf{C}^{\intercal}+\mathbf{C}\mathbf{C}^{\intercal}\mathbf{C}+\mathbf{C}(\mathbf{C}^{\intercal})^{2}$
$\displaystyle+\mathbf{C}^{\intercal}\mathbf{C}^{2}+\mathbf{C}^{\intercal}\mathbf{C}\mathbf{C}^{\intercal}+(\mathbf{C}^{\intercal})^{2}\mathbf{C}+(\mathbf{C}^{\intercal})^{3}\Big{)}$
$\displaystyle+\frac{\varsigma^{-2}}{8}\Big{(}2\mathbf{C}^{3}+3\mathbf{C}^{2}\mathbf{C}^{\intercal}+2\mathbf{C}\mathbf{C}^{\intercal}\mathbf{C}+3\mathbf{C}(\mathbf{C}^{\intercal})^{2}$
$\displaystyle+\mathbf{C}^{\intercal}\mathbf{C}^{2}+2\mathbf{C}^{\intercal}\mathbf{C}\mathbf{C}^{\intercal}+(\mathbf{C}^{\intercal})^{2}\mathbf{C}+2(\mathbf{C}^{\intercal})^{3}\Big{)}+\mathcal{O}(\mathbf{C}^{4})$
$\displaystyle=$
$\displaystyle\varsigma^{-2}\mathbf{I}-\frac{\varsigma^{-2}}{2}(\mathbf{C}+\mathbf{C}^{\intercal})+\frac{\varsigma^{-2}}{4}\left(\mathbf{C}^{\intercal}\mathbf{C}-\mathbf{C}\mathbf{C}^{\intercal}\right)$
$\displaystyle+\frac{\varsigma^{-2}}{8}\left(\mathbf{C}\left(\mathbf{C}^{\intercal}\mathbf{C}-\mathbf{C}\mathbf{C}^{\intercal}\right)+\left(\mathbf{C}^{\intercal}\mathbf{C}-\mathbf{C}\mathbf{C}^{\intercal}\right)\mathbf{C}^{\intercal})\right)+\mathcal{O}(\mathbf{C}^{4}).$
(77)

The submatrix of the Hessian for $\mathbf{y},\mathbf{x}$ couplings is

$\displaystyle\mathbf{H}_{yx}=$
$\displaystyle-\frac{\varsigma^{-2}}{4}\Big{(}\mathbf{C}_{ys}\mathbf{C}_{xs}^{\intercal}+\mathbf{C}_{ya}\mathbf{C}_{xa}^{\intercal}\Big{)}$
$\displaystyle+\frac{\varsigma^{-2}}{8}\Big{(}\mathbf{C}_{ya}(-\mathbf{C}_{aa}\mathbf{C}_{xa}^{\intercal}-\mathbf{C}_{as}\mathbf{C}_{xs}^{\intercal}-\mathbf{C}_{ax}\mathbf{C}_{xx}^{\intercal}+\mathbf{C}_{aa}^{\intercal}\mathbf{C}_{ax}-\mathbf{C}_{aa}^{\intercal}\mathbf{C}_{xa}^{\intercal}$
$\displaystyle-\mathbf{C}_{sa}^{\intercal}\mathbf{C}_{xs}^{\intercal}+\mathbf{C}_{xa}^{\intercal}\mathbf{C}_{xx}-\mathbf{C}_{xa}^{\intercal}\mathbf{C}_{xx}^{\intercal})$
$\displaystyle+\mathbf{C}_{ys}(-\mathbf{C}_{sa}\mathbf{C}_{xa}^{\intercal}-\mathbf{C}_{ss}\mathbf{C}_{xs}^{\intercal}+\mathbf{C}_{as}^{\intercal}\mathbf{C}_{ax}-\mathbf{C}_{as}^{\intercal}\mathbf{C}_{xa}^{\intercal}-\mathbf{C}_{ss}^{\intercal}\mathbf{C}_{xs}^{\intercal}+\mathbf{C}_{xs}^{\intercal}$
$\displaystyle\mathbf{C}_{xx}-\mathbf{C}_{xs}^{\intercal}\mathbf{C}_{xx}^{\intercal})+\mathbf{C}_{yy}(-\mathbf{C}_{ya}\mathbf{C}_{xa}^{\intercal}-\mathbf{C}_{ys}\mathbf{C}_{xs}^{\intercal}-\mathbf{C}_{sy}^{\intercal}\mathbf{C}_{xs}^{\intercal})$
$\displaystyle+\mathbf{C}_{sy}^{\intercal}(\mathbf{C}_{sa}\mathbf{C}_{xa}^{\intercal}+\mathbf{C}_{ss}\mathbf{C}_{xs}^{\intercal})+\mathbf{C}_{yy}^{\intercal}(\mathbf{C}_{ya}\mathbf{C}_{xa}^{\intercal}+\mathbf{C}_{ys}\mathbf{C}_{xs}^{\intercal})\Big{)}+\mathcal{O}(\mathbf{C}^{4})$
(78)

###  B.2 Solenoidal flows

We can describe the surprise of the system and its gradient as

$\displaystyle\textfrak{J}(\mathbf{z})=$ $\displaystyle-\log p(\mathbf{z})$
$\displaystyle=$
$\displaystyle\frac{1}{2}(\mathbf{z}-\mathbf{m}_{z})^{\intercal}\boldsymbol{\Sigma}_{zz}^{-1}(\mathbf{z}-\mathbf{m}_{z})+\frac{k}{2}\log\left(2\pi\right)+\frac{1}{2}\log|\boldsymbol{\Sigma}_{zz}|,$
(79) $\displaystyle\nabla_{\mathbf{z}}\textfrak{J}(\mathbf{z})=$
$\displaystyle\left(\boldsymbol{\Sigma}^{*}\right)^{-1}(\mathbf{z}-\mathbf{m}_{z}).$
(80)

where we know from Eq. 25 that
$\boldsymbol{\Sigma}^{*}=-2\mathbf{J}^{-1}\boldsymbol{\Gamma}-\mathbf{J}^{-1}\boldsymbol{\Sigma}^{*}\mathbf{J}^{\intercal}$.

At the NESS, the solution can be described in two terms consistent with a SDE
decomposition:

$\displaystyle\mathbf{J}(\mathbf{z}-\boldsymbol{\rho})=$
$\displaystyle(\mathbf{Q}-\boldsymbol{\Gamma})(-2\mathbf{J}^{-1}\boldsymbol{\Gamma}-\mathbf{J}^{-1}\boldsymbol{\Sigma}^{*}\mathbf{J}^{\intercal})^{-1}(\mathbf{z}-\boldsymbol{\rho}),$
(81) $\displaystyle\mathbf{Q}=$
$\displaystyle\boldsymbol{\Gamma}-\mathbf{J}(2\mathbf{J}^{-1}\boldsymbol{\Gamma}+\mathbf{J}^{-1}\boldsymbol{\Sigma}^{*}\mathbf{J}^{\intercal})=-\boldsymbol{\Sigma}^{*}\mathbf{J}^{\intercal}-\boldsymbol{\Gamma}=\boldsymbol{\Gamma}+\mathbf{J}\boldsymbol{\Sigma}^{*}.$
(82)

Rearranging the terms in the equation above to substitute
$\boldsymbol{\Sigma}^{*}=\mathbf{J}^{-1}(\mathbf{Q}-\boldsymbol{\Gamma})$ we
obtain the equivalence

$\displaystyle\mathbf{J}\mathbf{Q}+\mathbf{Q}\mathbf{J}^{\intercal}=\mathbf{J}\boldsymbol{\Gamma}-\boldsymbol{\Gamma}\mathbf{J}^{\intercal},$
(83)

which again takes the form of a continuous Lyapunov equation. This equation
can be solve numerically, or analytically using a power series expansion of
$\mathbf{J}=-\mathbf{I}+\mathbf{C}$, expanding the expression

$\displaystyle\mathbf{Q}=\frac{1}{2}(\mathbf{C}\boldsymbol{\Gamma}-\boldsymbol{\Gamma}\mathbf{C}^{\intercal})+\frac{1}{2}(\mathbf{C}\mathbf{Q}+\mathbf{Q}\mathbf{C}^{\intercal}),$
(84)

into the power series

$\displaystyle\mathbf{Q}$
$\displaystyle=\frac{1}{2}(\mathbf{C}\boldsymbol{\Gamma}-\boldsymbol{\Gamma}\mathbf{C}^{\intercal})+\frac{1}{4}\left(\mathbf{C}^{2}\boldsymbol{\Gamma}-\boldsymbol{\Gamma}(\mathbf{C}^{\intercal})^{2}\right)$
$\displaystyle+\frac{1}{8}\left(\mathbf{C}^{3}\boldsymbol{\Gamma}+\mathbf{C}^{2}\boldsymbol{\Gamma}\mathbf{C}^{\intercal}-\mathbf{C}\boldsymbol{\Gamma}(\mathbf{C}^{\intercal})^{2}-\boldsymbol{\Gamma}(\mathbf{C}^{\intercal})^{3}\right)+\mathcal{O}(\mathbf{C}^{4}).$
(85)

###  B.3 Mapping between the most likely internal and external states

If $\mathbf{z}$ can be divided into blanket states $\mathbf{b}$, and
internal/external states $\mathbf{h}=\\{\mathbf{x},\mathbf{y}\\}$, the
conditional distribution given blanket states is a Normal distribution with
moments

$\displaystyle\mathbf{m}_{h}(\mathbf{b}_{t})=$
$\displaystyle\boldsymbol{\rho}_{h}+\boldsymbol{\Sigma}^{*}_{hb}\left(\boldsymbol{\Sigma}^{*}_{bb}\right)^{g}(\mathbf{b}_{t}-\boldsymbol{\rho}_{b}),$
(86) $\displaystyle\boldsymbol{\Sigma}_{hh}(\mathbf{b}_{t})=$
$\displaystyle\boldsymbol{\Sigma}^{*}_{hh}-\boldsymbol{\Sigma}^{*}_{hb}\left(\boldsymbol{\Sigma}^{*}_{bb}\right)^{g}\boldsymbol{\Sigma}^{*}_{bh}.$
(87)

Note that the Markov blanket imposes that
$\boldsymbol{\Sigma}_{xy}(\mathbf{b}_{t})=\boldsymbol{\Sigma}_{yx}(\mathbf{b}_{t})=0$,
since $\mathbf{x},\mathbf{y}$ are conditionally independent.

Similarly, we can decompose $\mathbf{h}$ into an internal and external state

$\displaystyle\mathbf{m}_{x}(\mathbf{b}_{t})=$
$\displaystyle\boldsymbol{\rho}_{x}+\boldsymbol{\Sigma}^{*}_{xb}\left(\boldsymbol{\Sigma}^{*}_{bb}\right)^{g}(\mathbf{b}_{t}-\boldsymbol{\rho}_{b}),$
(88) $\displaystyle\mathbf{m}_{y}(\mathbf{b}_{t})=$
$\displaystyle\boldsymbol{\rho}_{y}+\boldsymbol{\Sigma}^{*}_{yb}\left(\boldsymbol{\Sigma}^{*}_{bb}\right)^{g}(\mathbf{b}_{t}-\boldsymbol{\rho}_{b}).$
(89)

If the mapping from $\mathbf{b}$ to $\mathbf{x}$ is injective (for this
$n_{x}\geq n_{b}$ is required), the first term can be rearranged into

$\displaystyle\mathbf{b}_{t}-\boldsymbol{\rho}_{b}=\boldsymbol{\Sigma}^{*}_{bb}\left(\boldsymbol{\Sigma}^{*}_{xb}\right)^{-1}\left(\mathbf{m}_{x}(\mathbf{b}_{t})-\boldsymbol{\rho}_{x}\right),$
(90)

yielding

$\displaystyle\mathbf{m}_{y}(\mathbf{b}_{t})=$
$\displaystyle\mathbf{m}_{y}+\boldsymbol{\Sigma}^{*}_{yb}\left(\boldsymbol{\Sigma}^{*}_{bb}\right)^{g}\boldsymbol{\Sigma}^{*}_{bb}\left(\boldsymbol{\Sigma}^{*}_{xb}\right)^{-1}\left(\mathbf{m}_{x}(\mathbf{b}_{t})-\boldsymbol{\rho}_{x}\right),$
$\displaystyle=$
$\displaystyle\boldsymbol{\rho}_{y}+\boldsymbol{\Sigma}^{*}_{yb}\left(\boldsymbol{\Sigma}^{*}_{xb}\right)^{-1}\left(\mathbf{m}_{x}(\mathbf{b}_{t})-\boldsymbol{\rho}_{x}\right)\equiv\sigma(\mathbf{m}_{x}(\mathbf{b}_{t})),$
(91)

which is a linear mapping of $\mathbf{m}_{x}(\mathbf{b}_{t})$ to
$\mathbf{m}_{y}(\mathbf{b}_{t})$.

$\displaystyle\mathbf{m}_{x}(\mathbf{b}_{t})=$
$\displaystyle\mathbf{m}_{x}+\boldsymbol{\Sigma}^{*}_{xb}\left(\boldsymbol{\Sigma}^{*}_{bb}\right)^{g}\boldsymbol{\Sigma}^{*}_{bb}\left(\boldsymbol{\Sigma}^{*}_{yb}\right)^{-1}\left(\mathbf{m}_{y}(\mathbf{b}_{t})-\mathbf{m}_{y}(\mathbf{b}_{t})\right)$
$\displaystyle=$
$\displaystyle\mathbf{m}_{x}+\boldsymbol{\Sigma}^{*}_{xb}\left(\boldsymbol{\Sigma}^{*}_{yb}\right)^{-1}\left(\mathbf{m}_{y}(\mathbf{b}_{t})-\mathbf{m}_{y}(\mathbf{b}_{t})\right)\equiv\sigma^{-1}(\mathbf{m}_{y}(\mathbf{b}_{t})).$
(92)

In the NESS limit

$\displaystyle\frac{\mathrm{d}\mathbf{m}_{y}(\mathbf{b}_{t})}{\mathrm{d}t}=$
$\displaystyle\boldsymbol{\Sigma}^{*}_{yb}\left(\boldsymbol{\Sigma}^{*}_{xb}\right)^{-1}\frac{\mathrm{d}\mathbf{m}_{x}(\mathbf{b}_{t})}{\mathrm{d}t}=\nabla_{\mathbf{m}_{x}}\sigma(\mathbf{m}_{x}(\mathbf{b}_{t}))\frac{\mathrm{d}\mathbf{m}_{x}(\mathbf{b}_{t})}{\mathrm{d}t}.$
(93)

Similarly, if the mapping from $\mathbf{b}$ to $\mathbf{y}$ is injective (for
this $n_{y}\geq n_{b}$ is required), we can invert the relation

$\displaystyle\frac{\mathrm{d}\mathbf{m}_{x}(\mathbf{b}_{t})}{\mathrm{d}t}=$
$\displaystyle\boldsymbol{\Sigma}^{*}_{xb}\left(\boldsymbol{\Sigma}^{*}_{yb}\right)^{-1}\frac{\mathrm{d}\mathbf{m}_{y}(\mathbf{b}_{t})}{\mathrm{d}t}=\nabla_{\mathbf{m}_{y}}\sigma^{-1}(\mathbf{m}_{y}(\mathbf{b}_{t}))\frac{\mathrm{d}\mathbf{m}_{y}(\mathbf{b}_{t})}{\mathrm{d}t}.$
(94)

##  Appendix C How informative are conditional marginal flows? A minimal
example

In order to illustrate the debate about how informative are average flows, we
study a simple two-dimensional linear stochastic system

$\displaystyle\frac{\mathrm{d}y_{t}}{\mathrm{d}t}=$ $\displaystyle
J_{yy}y_{t}+J_{yb}b_{t}+\omega_{y,t},$ (95)
$\displaystyle\frac{\mathrm{d}b_{t}}{\mathrm{d}t}=$ $\displaystyle
J_{by}y_{t}+J_{bb}b_{t}+\omega_{b,t},$ (96)

which can be studied as the Langevin dynamics described in B. We select a
variance of $\omega_{y,t},\omega_{b,t}$, $\varsigma^{2}=1$, and

$\displaystyle\mathbf{J}=\begin{bmatrix}-1&1.2\\\ -1&-1\end{bmatrix}.$ (97)

The resulting covariance matrix is

$\displaystyle\boldsymbol{\Sigma}=\begin{bmatrix}0.9545&0.04545\\\
0.04545&1.05455\end{bmatrix}.$ (98)

Fig. 1 displays an example of a trajectory of this system.

###  C.1 Conditional marginal flows capture partial tendencies, not real
behaviour

Flows in the system are described as

$\displaystyle f_{y}(y_{t},b_{t})=$ $\displaystyle J_{yy}y_{t}+J_{yb}b_{t},$
(99) $\displaystyle f_{b}(y_{t},b_{t})=$ $\displaystyle
J_{by}y_{t}+J_{bb}b_{t}.$ (100)

Following the results in B, conditional average flows can be described as

$\displaystyle\left\langle f_{y}(y_{t},b_{t})\right\rangle_{b_{t}}=$
$\displaystyle J_{yy}m_{y}(b_{t})+J_{yb}b_{t}$ $\displaystyle=$
$\displaystyle\left(J_{yy}\Sigma^{*}_{yb}/\Sigma^{*}_{bb}+J_{yb}\right)b_{t},$
(101) $\displaystyle\left\langle f_{y}(y_{t},b_{t})\right\rangle_{y_{t}}=$
$\displaystyle\left(J_{yy}+J_{yb}\Sigma^{*}_{by}/\Sigma^{*}_{yy}\right)y_{t},$
(102) $\displaystyle\left\langle f_{b}(y_{t},b_{t})\right\rangle_{b_{t}}=$
$\displaystyle\left(J_{by}\Sigma^{*}_{yb}/\Sigma^{*}_{bb}+J_{bb}\right)b_{t},$
(103) $\displaystyle\left\langle f_{b}(y_{t},b_{t})\right\rangle_{y_{t}}=$
$\displaystyle\left(J_{by}+J_{bb}\Sigma^{*}_{by}/\Sigma^{*}_{yy}\right)y_{t}.$
(104)

In Fig. 2.A we describe the flow structure in the system, which describes a
spiral behaviour due to the non-equilibrium tendencies in the system. This is
an example of the solenoidal flows captured by the matrix $\mathbf{Q}$.

In contrast, Fig. 2.B captures the marginal flow structure when variable $b$
is fixed. This separates tendencies in the system, resulting in a diverging
tendency for variable $y$ (diverging to $-\infty$ for negative $b$ and to
$+\infty$ for positive $b$. Similarly, combining cross marginal flows (fix $b$
for the marginal flow of $y$ and vice versa) result in another combination of
partial tendencies as shown in Fig. 2.C. In this case, the rotation of
solenoidal flows is captured, but no the attractor that structures this
rotation into an spiral behaviour.

In sum, conditonal marginal flows do not capture the behaviour of the system.
Even if the flows point in some direction (e.g. free energy minimization) we
cannot conclude that this is equivalent to the system behaving ‘as if’
following that particular direction. An alternative illustration of this
problem can be described by writing the conditional average flow of $y$ in
terms of the most likely state $m_{y}(b)$

$\displaystyle\left\langle f_{y}(y_{t},b_{t})\right\rangle_{b_{t}}=$
$\displaystyle\left(J_{yy}+J_{yb}\Sigma^{*}_{by}/\Sigma^{*}_{yy}\right)m_{y}(b_{t}).$
(105)

Where the term $\left(J_{yy}+J_{yb}\Sigma^{*}_{by}/\Sigma^{*}_{yy}\right)$ is
a positive constant (see Fig. 3.B). Following the logic of a relaxed
interpretation of the FEP, this could be interpreted as a gradient ascent on
$m_{y}^{2}(b_{t})$. However, as we see in the Fig. 2.A the real behaviour of
the system will display a global attractor at $y=b=0$ with solenoidal
couplings, eventually minimizing $m_{y}^{2}(b_{t})$.

Figure 1: Evolution of a bivariate linear stochastic system. Example trajectory of the system described by Eq. 96 with parameters $\varsigma=1$, $J_{yy}=J_{bb}=-1$, $J_{yb}=1.2$, $J_{by}=-1$.  A | B | C
---|---|---
|  |
Figure 2: Flow structure vs conditional average flow structure of the system.
Comparison between the flow structure in the system displaying a spiral
behaviour (A), the conditional average flows when variable $b$ is fixed,
displaying divergent tendencies in $y$ and convergent tendencies in $b$ (B),
and the cross conditional average flow when variable $b$ is fixed for $f_{y}$
and variable $y$ is fixed for $f_{b}$, displaying rotation without an
attractor (B).

###  C.2 Rates of conditional averages are different to conditional marginal
flows

The second issue we illustrate in this simple system is the difference between
conditional marginal flows and the dynamics of the most likely states.

From B we can derive

$\displaystyle\frac{\mathrm{d}m_{y}(b_{t})}{\mathrm{d}t}=$
$\displaystyle{\Sigma}^{*}_{yb}/{\Sigma}^{*}_{bb}\frac{\mathrm{d}b_{t}}{\mathrm{d}t}={\Sigma}^{*}_{yb}/{\Sigma}^{*}_{bb}\left(J_{by}y_{t}+J_{bb}b_{t}+\varsigma\omega_{b,t}\right).$
(106)

For a fixed $b$, this variable is distributed as a Normal distribution
$\mathcal{N}(u,s)$, with mean and standard deviation

$\displaystyle u=$
$\displaystyle{\Sigma}^{*}_{yb}/{\Sigma}^{*}_{bb}J_{by}m_{y}(b),$ (107)
$\displaystyle s=$
$\displaystyle{\Sigma}^{*}_{yb}/{\Sigma}^{*}_{bb}\left(J_{by}(\Sigma^{*}_{yy}-\Sigma^{*}_{yb}\Sigma^{*}_{by}/\Sigma^{*}_{bb})+\sqrt{2}\varsigma
s_{\omega}\right),$ (108)

where $s_{\omega}$ is the standard deviation of the noise introduced in the
Langevin dynamics (described as $\sqrt{\frac{1}{\mathrm{d}t}}$).

For a range of values of $b$, Fig. 3 captures the conditional average flows
$\left\langle f_{y}(y_{t},b_{t})\right\rangle_{b_{t}}$ (solid dark line), and
the distribution of the derivatives of the conditional average state
$\frac{\mathrm{d}m_{y}(b_{t})}{\mathrm{d}t}$ (the mean is represented by the
dashed line and error bars of three standard deviations by the light area). As
we can observe, the dependency with $m_{y}(b)$ is positive with respect to
$\left\langle f_{y}(y_{t},b_{t})\right\rangle_{b_{t}}$ (indicating a gradient
ascent on $m_{y}^{2}(b)$) while the true dynamics is captured by the negative
dependency between $\frac{\mathrm{d}m_{y}(b_{t})}{\mathrm{d}t}$ and $m_{y}(b)$
(indicating a gradient descent on $m_{y}^{2}(b)$). Note that for different
parameters the sign of the slope of $\left\langle
f_{y}(y_{t},b_{t})\right\rangle_{b_{t}}$ can change (being negative or
positive) but the slope of $\frac{\mathrm{d}m_{y}(b_{t})}{\mathrm{d}t}$ is
always negative, given the presence of a global attractor. This shows how,
even in very simple examples, these quantities can have radically different
behaviours and that the conditional average flow does not necessarily capture
the true behaviour of a system.

A | B
---|---
|
Figure 3: Rate of conditional averages versus conditional average flows. The
conditional average flow of $y$ fixing $b$, $\left\langle
f_{y}(y_{t},b_{t})\right\rangle_{y_{t}}$, is represented by the solid dark
line. The rate of change of the conditional average
$\frac{m_{y}(b_{t})}{\mathrm{d}t}$ is represented by the dashed line
(representing its average value for a given $b$, and the light area
(representing error bars for 3 standard deviations for a fixed $b$.). These
quantities are represented with respect to a $b$ axis (A) and a $m_{y}(b)$
axis (B), as one is a linear transformation of the other.