Bayesian mechanics of perceptual inference and motor control in the brain

The free energy principle (FEP) in the neurosciences stipulates that all viable agents induce and minimize informational free energy in the brain to fit their environmental niche. In this study, we continue our effort to make the FEP a more physically principled formalism by implementing free energy minimization based on the principle of least action. We build a Bayesian mechanics (BM) by casting the formulation reported in the earlier publication (Kim in Neural Comput 30:2616–2659, 2018, 10.1162/neco_a_01115) to considering active inference beyond passive perception. The BM is a neural implementation of variational Bayes under the FEP in continuous time. The resulting BM is provided as an effective Hamilton’s equation of motion and subject to the control signal arising from the brain’s prediction errors at the proprioceptive level. To demonstrate the utility of our approach, we adopt a simple agent-based model and present a concrete numerical illustration of the brain performing recognition dynamics by integrating BM in neural phase space. Furthermore, we recapitulate the major theoretical architectures in the FEP by comparing our approach with the common state-space formulations.


Introduction
The free energy principle (FEP) in the field of neurosciences rationalizes that all viable organisms cognize and behave in the natural world by calling forth the probabilistic models in their neural system-the brain-in a manner that ensures their adaptive fitness (Friston 2010a). The neurobiological mechanism that endows an organism's brain-the neural observer-with this ability is theoretically framed into an inequality that weighs two information-theoretical measures: surprisal and informational free energy (IFE) (see, for a review, Buckley et al. 2017). The surprisal provides a measure of the atypicality of an environmental niche, and the IFE is the upper bound of the surprisal. The inequality enables a cognitive agent to minimize the IFE as a variational objective Communicated by Karl Friston. B Chang Sub Kim cskim@jnu.ac.kr 1 Department of Physics, Chonnam National University, Gwangju 61186, Republic of Korea function indirectly instead of the intractable surprisal. 1 The minimization corresponds to inferring the external causes of afferent sensory data, which are encoded as a probability density at the sensory interface, e.g., sensory organs. The brain of an organism neurophysically performs the Bayesian computation of minimizing the induced variational IFE; this is termed as recognition dynamics (RD), which emulates, under the Laplace approximation (Friston et al. 2007), the predictive coding scheme of message processing or recognition (Rao and Ballard 1999;Bogacz 2017). The neuronal self-organization in vitro under the FEP was studied recently at the level of single neuron responses (Isomura et al. 2015;Isomura and Friston 2018). Owing to its explanatory power of perception, learning, and behavior of living organisms within a framework, it is suggested a promising unified biological principle (Friston 2010a(Friston , 2013Colombo and Wright 2018;Ramstead et al. 2018).
The neurophysical mechanisms of the abductive inference in the brain are yet to be understood; therefore, researchers mostly rely on information-theoretic concepts (Elfwing et al. 2016;Ramstead et al. 2019;Kuzma 2019;Shimazaki 2019;Kiefer 2020, Sanders et al. 2020. The FEP facilitates dynamic causal models in the brain's generalized-state space (Friston 2008b;Friston et al. 2010b), which pose a mixed discrete-continuous Bayesian filtering (Jazwinski 1970;Balaji and Friston 2011). In this work, we consider that the brain confronts the continuous influx of stochastic sensations and conducts the Bayesian inversion of inferring external causes in the continuous state representations. Biological phenomena are naturally continuous spatiotemporal events; accordingly, we suggest that the continuous-state approaches used to describe cognition and behavior are better suited than discrete-state descriptions for studying the perceptual computation in the brain.
Recently, we carefully evaluated the FEP while clarifying technical assumptions that underlie the continuous statespace formulation of the FEP (Buckley et al. 2017). A full account of the discrete-state formulation complementary to our formulation can be found in (Da Costa et al. 2020a). In a subsequent paper (Kim 2018), we reported a different variational scheme that the Bayesian brain may utilize in conducting inference. In particular, by postulating that "surprisal" plays the role of a Lagrangian in theoretical mechanics (Landau and Lifshitz 1976;Sengupta et al. 2016), we worked a plausible computational implementation of the FEP by utilizing the principle of least action. We believed that although the FEP relies on Bayesian abductive computation, it must be properly formulated conforming to the physical principles and laws governing the matter comprising the brain. To this end, we proposed that any process theory of the FEP ought to be based on the full implication of the inequality (Kim 2018) where ϕ and ϑ collectively denote the sensory inputs and environmental hidden states, respectively. The integrand on the left-hand side (LHS) of the preceding equation − ln p(ϕ) is the aforementioned surprisal, which measures the "selfinformation" contained in the sensory density p(ϕ) (Cover and Thomas 2006), and F on the right-hand side (RHS) is the variational IFE defined as which encapsulates the recognition (R-) density q(ϑ) and the generative (G-) density p(ϑ, ϕ) (Buckley et al. 2017). While the G-density represents the brain's belief (or assumption) of sensory generation and hidden environmental dynamics, the R-density is the brain's current estimate of the environmental cause of the sensory perturbation. The G-and R-densities together induce variational IFE when receptors at the brainenvironment interface are excited by sensory perturbations. According to Eq. (1), the FEP articulates that the brain minimizes the upper bound of the sensory uncertainty, which is a long-term surprisal. We identify this bound as an informational action (IA) within the scope of the mechanical principle of least action (Landau and Lifshitz 1976). Then, by complying with the revised FEP, we formulate the Bayesian mechanics (BM) that executes the RD in the brain. The RD neurophysically performs the computation for minimizing the IA when the neural observer encounters continuous streams of sensory data. The advantage of our formulation is that the brain and the environmental states are specified using only bare continuous variables and their first-order derivatives (velocities or equivalent momenta). The momentum variables represent prediction errors, which quantify the discrepancy between an observed input and its top-down belief of a cognitive agent in the predictive coding language (Huang and Rao 2011;de Gardelle et al. 2013;Kozunov et al. 2020).
The goal of this work is to cast our previous study to include the agent's motor control, which acts on the environment to alter sensory inputs. 2 Previously, by utilizing the principle of least action, we focused on the formulation of perceptual dynamics for the passive inference of static sensory inputs (Kim 2018) without incorporating motor control for the active perception of nonstationary sensory streams. Here, we apply our approach to the problem of active inference derived from the FEP (Friston et al. 2009(Friston et al. , 2010c(Friston et al. , 2011a, which proposes that organisms can minimize the IFE by altering sensory observations when the outcome of perceptual inference alone is not in accordance with the internal representation of the environment (Buckley et al. 2017). Living systems are endowed with the ability to adjust their sensations via proprioceptive feedback, which is attributed to an inherited trait of all motile animals embodied in the reflex pathways (Tuthill and Azim 2018). In this respect, motor control is considered an inference of the causes of motor signals encoded as prediction errors at proprioceptors, and motor inference is realized at the spinal level by classical reflex arcs (Friston 2011b;Adams et al. 2013). Our formulation evinces time-dependent driving terms in the obtained BM, which arise from sensory prediction errors, as control (motor) signals. Accordingly, the BM bears resemblance to the deterministic control derived from Pontryagin's maxi-mum principle in optimal control theory (Todorov 2007). In this work, we consider the agent's locomotion for action inference only implicitly: our formulation focuses on the implementation of the control signal (or commands) at the neural level of description and not at the behavioral level of biological locomotion; accordingly, the additional minimization mechanism of the IA inferring optimal control was not explicitly handled, which is left as a future work. There are other systematic approaches that try to relate the active inference formalism to the existing control theories (Baltieri and Buckley 2019;Millidge et al. 2020a;Da Costa et al. 2020c).
Technically, a variation in the IA yields the BM that computes Bayesian inversion, which is given as a set of coupled differential equations for the brain variables and their conjugate momenta. The brain variables are ascribed to the brain's representation of the environmental states, and their conjugate momenta are the combined prediction errors of the sensory data and the rate of the state representations. The neural computation of active inference corresponds to the BM integration and is subject to nonautonomous motor signals. The obtained solution results in optimal trajectories in the perceptual phase space, which yields a minimum accumulation of the IFE over continuous time, i.e., a stationary value of the IA. Our IA is identical to that of "free action" defined in the Bayesian filtering schemes (Friston et al. 2008c). When the minimization of free action is formulated in the generalized filtering scheme (Friston et al. 2010b), two approaches are akin to each other such that both assume the Laplace-encoded IFE as a mechanical Lagrangian. The difference lies in their mathematical realization of minimization: our approach applies the principle of least action in classical mechanics, while generalized filtering uses the gradient descent method in the generalized state space, where generalized states are interpreted as a solenoidal gradient flow in a nonequilibrium steady state.
The remainder of this paper is organized as follows. In Sect. 2, we unravel some of the theoretical details in the formulation of the FEP. In Sect. 3, we formulate the BM of the sensorimotor cycle by utilizing the principle of least action. Then, in Sect. 4, we present a parsimonious model for a concrete manifestation of our formulation. Finally, in Sect. 5, we provide the concluding remarks.

Recapitulation of technical developments
Here, we recapitulate theoretical architectures in the continuousstate formulation under the FEP while discussing technical features that distinguish our formulation from prevailing state-space approaches.

Perspective on generalized states
The Bayesian filtering formalism of the FEP adopts the concept of the generalized motion of a dynamical object by defining its mechanical state beyond position and velocity (momentum). The generalized states of motion are generated by recursively taking time derivatives of the bare states. A point in the hyperspace defined by the generalized states is interpreted as an instantaneous trajectory. This notion provides an essential theoretical basis for ensuring an equilibrium solution of the RD in the conventional formulation of the FEP (Kim 2018); it is commonly employed by researchers (Parr and Friston 2018;Baltieri and Buckley 2019).
The motivation behind the generalized coordinates of motion is to describe noise correlation in the generative processes beyond white noise (Wiener process), and thus, to provide a more detailed specification of the dynamical states (Friston 2008a, b;Friston et al. 2010b). The mathematical theory of a quasi-Markovian process undergirds this formulation, which describes general stochastic dynamics with a colored-noise correlation with a finite-dimensional Markovian equation in an extended state space by adding auxiliary variables (Pavliotis 2014). State-space augmentation in terms of generalized coordinates may be considered a special realization of the Pavliotis formalism. The stateextension procedure adopts some specific approximations, such as the local linearization procedure developed in nonlinear time-series analysis (Ozaki 1992).
From the physics perspective, higher-order states possess a different dynamical status in comparison with Newtonian mechanical states, specified only by position (bare order) and velocity (first order). A change in the Newtonian states is caused by a force that specifies acceleration (second order) (Landau and Lifshitz 1976). Although there are no "generalized forces" causing the jerk (third order), snap (fourth order), etc., the jerk can be measured phenomenologically by observing a change in acceleration. This induces all higherorder states to the kinematic level. Another perspective is whether update equations in terms of generalized coordinates are equivalent to the Pavliotis' quasi-Markovian description. Auxiliary variables in Pavliotis' analysis are not generated by the recursive temporal derivatives of a bare state. The generalized phase space considered in (Kerr and Graham 2000) is also spanned in terms of canonical displacement and momentum variables. A further in-depth analysis is required.
Our formulation does not employ the generalized states, but instead, it follows the normative rules in specifying generative models (Kim 2018). The derived BM performs the brain's Bayesian inference in terms of only the bare brain variable and its conjugate momentum in phase space, and not in an extended state space. Accordingly, our formulation is restricted to the white noise in the generative processes; however, it provides a natural approach to determine the equi-librium solutions of the BM (see Sect. 4). For the general brain models described by many brain variables, the brain's BM can be set up in multi-dimensional phase space, which is distinctive from the state-space augmentation in the generalized coordinate formulation (see Sect. 2.3).

Continuous state implementation of recognition dynamics (RD)
The conventional FEP employs the gradient-descent minimization of the variational IFE by the brain's internal states. To incorporate the time-varying feature of sensory inputs, the method distinguishes the path of a mode and the mode of a path in the generalized state space (Friston 2008b;Friston et al. 2008cFriston et al. , 2010b. This theoretical construct intuitively considers the nonequilibrium dynamics of generalized brain states as drift-diffusion flows that locally conserve the ensemble density in the hyperspace of the generalized states (Friston et al. 2012b;Friston 2019). Mathematically, the gradient descent formulation is based on the general idea for a fast and efficient convergence and it ensures that formulations reach a sophisticate level by incorporating the Riemannian metric in information geometry (Amari 1998;Surace et al. 2020); the idea is applied to the FEP (Sengupta and Friston 2017;Da Costa et al. 2020b).
In our proposed formulation, we replace the gradient descent scheme with the standard mechanical formulation of the least action principle (Kim 2018). However, there is a disadvantage in that we incorporate only the Gaussian white noise in the generative processes of the sensory data and environmental dynamics [see Sect. 2.3]. The resulting novel RD described by an effective Hamiltonian mechanics entails optimal trajectories but no single fixed points in the canonical state (phase) space, which provides an estimate of the minimum sensory uncertainty, i.e., the average surprisal over a finite temporal horizon. The phase space comprises the positions (predictions) and momenta (prediction errors) of the brain's representations of the causal environment.
Our implementation of the minimization procedure is an alternative to the gradient descent algorithms in the FEP. A crucial difference between the two approaches is that while the gradient descent scheme searches for an instantaneous trajectory representing a local minimum on the IFE landscape in the multidimensional generalized state space, our theory determines an optimal trajectory minimizing the continuoustime integral of the IFE in two-dimensional phase space for a single variable problem.

Treatment of noise correlations
The FEP requires the brain's internal model of the G-density p(ϕ, ϑ) encapsulating the likelihood p(ϕ|ϑ) and prior p(ϑ). The likelihood density is determined by the random fluctu-ation in the expected sensory-data generation, and the prior density is determined by that in the believed environmental dynamics. The brain encounters sensory signals on a timescale, which is often shorter than the correlation time of the random processes (Friston 2008a); accordingly, in general, the noises embrace a non-Markovian stochastic process with an intrinsic temporal correlation that surmounts the ideal white-noise stochasticity. Conventional formulations (Friston 2008b; Friston et al. 2010b) consider that colored noises are analytic (i.e., differentiable) to allow correlation between the distinct dynamical orders of the continuous states. In practice, to furnish a closed dynamics for a finite number of variables, the recursive equations of motion for the continued generalized states need to be truncated at an arbitrary embedding order.
Our formulation considers the BM in the brain in terms of the standard Newtonian (Hamiltonian) construct; the drawback is that our theory does not explore the nature of temporal correlation in the assumed Gaussian noises in the generative processes. Accordingly, our generative models assume and account for the white noise describing the Wiener processes. The delta-correlated white noise is mathematically singular; they need to be smoothed to describe fast biophysical processes. There are approaches in stochastic theories that formulate non-Markovian processes with colored noises without resorting to generalized states of motion (van Kampen 1981;Fox 1987;Risken 1989;Moon and Wettlaufer 2014), which are not discussed here.
Instead, we discuss an approach to extend the phase-space dimension for the white noise processes. At the level of the Hodgkin-Huxley description of the biophysical brain dynamics, the membrane potential, gating variables, and ionic concentrations are relevant coarse-grained brain variables (Hille 2001). Thus, if one employs the fluctuating Hodgkin-Huxley models with Gaussian white noises as neurophysically plausible generative models (Kim 2018), one can proceed with our Lagrangian (equivalently, Hamiltonian) approach to formulate the RD in an extended phase space. Such a state-space augmentation is different from and alternative to that in terms of the generalized coordinates of motion [see Sect. 2.1], while accommodating only deltacorrelated noises.

Lagrangian formulation of Bayesian mechanics (BM)
The (classical) "action" is defined as an ordinary timeintegral of the Lagrangian for an arbitrary trajectory (Landau and Lifshitz 1976). Our formulation of the BM proposes the Laplace-encoded IFE-an upper bound on the sensory surprisal-as an informational Lagrangian and hypothesizes the time-integral of the IFE-an upper bound on the sensory Shannon uncertainty-as an informational action (IA). By applying the principle of least action, we minimize the IA to find a tight bound for the sensory uncertainty and derive the BM that performs the brain's Bayesian inference of the external cause of sensory data. In turn, we cast the working BM in our formulation as effective Hamilton's equations of motion in terms of position and momentum in phase space. Meanwhile, the BM described in (Friston 2019) intuitively adopts the idea of Feynman's path integral formulation (Feynman and Hibbs 2005). The Feynman's path integral formulation extends the idea of classical action to quantum dynamics and provides an approach to determine the "propagator" that specifies the transition probability between initial and final states. The propagator is defined as a functional integral of the exponentiated action, which summates all possible trajectories connecting initial and final states. The description provided in (Friston 2019) identifies the propagator using the probability density over neural states, and it makes the connection to the Bayesian FEP. In this manner, the surprisal may be identified as a negative log of the steady-state density in nonequilibrium ensemble dynamics (Parr et al. 2020), which is governed by a Fokker-Plank equation. The generalized Bayesian filtering scheme (Friston et al. 2008c(Friston et al. , 2010b provides a continuous-state formulation of minimizing the surprisal, and it delivers the BM in terms of the generalized coordinates of motion using the concept of gradient flow. In some technical details, the Lagrangian presented in (Friston 2019), which is the integrand in the classical action, encloses two terms. They are the quadratic term arising from the state equation and the term involving a state-derivative (divergence in three dimension) of the force, which appears in the Langevin-type state equation. The former term is included in our Lagrangian but with an additional quadratic term from the observation equation. In contrast, the latter is not present in our Lagrangian, which is known to arise from the Stratonovich convention (Seifert 2012; Cugliandolo and Lecomte 2017).

Closure of the sensorimotor loop in active inference
The conventional FEP facilitates gradient descent minimization for the mechanistic implementation of active inference, which makes the motor-control dynamics available in the brain's RD (Friston et al. 2009(Friston et al. , 2010c(Friston et al. , 2011a. The gradientdescent scheme is mathematically expressed aṡ where a denotes an agent's motor variable and F represents the Laplace-encoded IFE by the biophysical brain variables (Buckley et al. 2017). An agent's capability of subjecting sensory inputs to motor control is considered a functional dependence ϕ = ϕ(a) in the environmental generative processes (Friston et al. 2009). According to Eq. (3), an agent performs the minimization by effectuating the sensory data dϕ/da and obtains the best result for motor inference wheṅ a = 0, where the condition ∂ F/∂ϕ = 0 must be met. Because ∂ F/∂ϕ produces terms proportional to the sensory prediction errors, the fulfillment of motor inference is equivalent to suppressing proprioceptive errors. Thus, motor control attempts to minimize prediction errors, while prediction errors convey motor signals for control dynamics; this forms a sensorimotor loop. Some subtle questions arise here regarding the dynamical status of the motor-control variable a: Eq. (3) evidently handles a as a dynamical state; however, the corresponding equation of motion governing its dynamics is not given in the environmental processes. Instead, the mechanism of motor control that vicariously alters the sensory-data generation is presumed (Friston et al. 2009). In addition, motor variables are represented as the active states of the brain, e.g., motor-neuron activities in the ventral horn of the spinal cord (Friston et al. 2010c); however, they are treated differently from other hidden-state representations.
Recall that the internal state representations are expressed as generalized states, whereas the active states are not.
In the following, we pose a semi-active inference problem that does not explicitly address optimal motor control (motor inference) in the RD but encompasses the motorcontrol signal as a time-dependent driving term arising from nonstationary prediction errors in the sensory-data cause.

Closed-loop dynamics of perception and motor control
The brain is not divided into sensory and motor systems. Instead, it is one inference machine that performs the closedloop dynamics of perception and motor control. Here, we develop a framework of active inference within the scope of the least action principle by employing the Laplace-encoded IFE as an informational Lagrangian. The environmental states ϑ undergo deterministic or stochastic dynamics by obeying physical laws and principles. Here, we do not explicitly consider their equations of motion because they are hidden from the brain's perspective, i.e., the brain as a neural observer does not possess direct epistemic access. Similarly, sensory data ϕ are physically generated by an externally hidden process at a sensory receptor, which constitutes the brain-environment interface. However, to emphasize the effect of an agent's motor control a on sensory generation, we facilitate the generative process of sensory data using an instantaneous mapping where h(ϑ, a) denotes the linear or nonlinear map of input generation and z gp represents the noise involved. Note that an agent's motor-control a is explicitly included in the generative map. However, the neural observer is not aware of how the sensory streams are effectuated by the agent's motion in the environment (Friston et al. 2010c). The FEP circumvents this epistemic difficulty by hypothesizing a formal homology between external physical processes and the corresponding internal models foreseen by the neural observer (Friston et al. 2010c). Upon receiving sensory-data influx, the brain launches R-density q(ϑ) to infer the external causes via variational Bayes. The R-density is the probabilistic representation of the environment, whose sufficient statistics are assumed to be encoded by neurophysical brain variables, e.g., neuronal activity or synaptic efficacy. When a fixed-form Gaussian density is considered for the R-density, which is called Laplace approximation, only the first-order sufficient statistic, i.e., the mean μ, is needed to specify the IFE effectively (Buckley et al. 2017). The brain continually updates the R-density using its internal dynamics, described here as a Langevin-type equation where f (μ) represents the brain's belief regarding the external dynamics encoded by a neurophysical driving mechanism of the brain variables μ, and w is random noise. The sensory perturbations at the receptors are predicted by the neural observer via the instantaneous mapping where the belief g(μ) is encoded by the internal variables, and z is the associated noise. Our sensory generative model provides a mechanism for sampling sensory data ϕ using the brain's active states a, which represent an external motor control embedded in Eq. (4). Note that Eq. (4) describes the environmental processes that generate sensory inputs ϕ, while its homolog "Eq. (6)" prescribes the brains' prior belief of ϕ that can be altered by the active states a. The instantaneous state of the brain μ, which is specified by Eq. (5), selects a particular R-density q(ϑ) when the brain seeks the true posterior (the goal of perceptual inference). The motor control fulfills the prior expectations by modifying the sensory generation via active-state effectuation at the proprioceptors. Through Laplace approximation (Buckley et al. 2017), the G-density p(ϕ, ϑ) is encoded in the brain as p = p(ϕ, μ), where the sensory stimuli ϕ are predicted by the neural observer μ via Eq. (6). Here, we argue that the physical sensory-recording process is conditionally independent of the brain's internal dynamics; however, the brain states must be neurophysically involved in computing the sensory prediction. In other words, from the physics perspective, the sensory perturbation ϕ at the interface is a source for exciting the neuronal activity μ. This observation renders the set of Eqs. (5) and (6) to be dynamically coupled, and not conditionally independent. We incorporate this conditional dependence into our formulation by introducing a statistical coupling via the covariance connection between the likelihood p(ϕ|μ) and prior p(μ) that together furnish the Laplace-encoded G-density.
For simplicity, we consider the stationary Gaussian processes for the bivariate variable Z as a column vector where w =μ − f (μ) and z = ϕ − g(μ), and we specify the Laplace-encoded G-density p(ϕ, μ) = p(ϕ|μ) p(μ) as where |Σ| and Σ −1 are the determinant and the inverse of the matrix Σ, respectively; Z T is the transpose of Z . The covariance matrix Σ for the above is given as where the stationary variances σ i (i = w, z) and the transient covariance φ are defined, respectively, as With the prescribed internal model of the brain for the G-density, the Laplace-encoded IFE can be specified as F(ϕ, μ) = − ln p(ϕ, μ) (for details, see Buckley et al. 2017). Then, it follows that where ρ denotes the correlation function defined as a normalized covariance Furthermore, we introduce notations m i (i = w, z) as which are precisions, scaled by the correlation in the conventional FEP. Next, as proposed in (Kim 2018), we identify F as an informational Lagrangian L within the scope of the principle of least action, and we define which is viewed as a function of μ andμ for the given sensory inputs ϕ(t), i.e., L = L(μ,μ; ϕ). Note that we dropped the last term in Eq. (8) when translating F into L because it can be expressed as a total time-derivative term that does not affect the resulting equations of motion (Landau and Lifshitz 1976). Then, the theoretical action S that effectuates the variational objective functional under the revised FEP is set up as The Euler-Lagrange equation of motion, which determines the trajectory μ = μ(t) for a given initial condition μ (0), is derived by minimizing the action δS ≡ 0. Equivalently, the equations of motion can be considered in terms of the position μ and its conjugate momentum p, instead of the position μ and velocityμ. We used the terms position and velocity as a metaphor to indicate dynamical variables μ andμ, respectively. For this, we need to convert Lagrangian L into Hamiltonian H by performing a Legendre transformation where p denotes the canonical momentum conjugate to μ, which is calculated from L as After some manipulation, the functional form of H can be obtained explicitly as where we indicated its dependence on the sensory influx ϕ.
In addition, the terms T and V on the RHS are defined as Here, T and V represent the kinetic and potential energies, respectively, which define the informational Hamiltonian of the brain. Similarly, m w and m z represent the neural inertial masse as a metaphor. Unlike that in standard mechanics, the second term in the expression for kinetic energy is dependent on linear momentum and position. We generate the Hamilton equations of motion, which are equivalent to the Lagrange equation usinġ As described below, Hamilton's equations are better suited for our purposes because they specify the RD as coupled first-order differential equations of the brain state μ and its conjugate momentum p. In contrast, the Lagrange equation is a second-order differential equation of the state variable (Landau and Lifshitz 1976). The results arė where parameters α, β, and γ have been, respectively, defined for notational convenience as where κ denotes the tuning parameter to spawn stability. In Eqs. (17) and (18), we defined the notation ϕ as It measures the discrepancy between the adjustable sensory input ϕ by an agent's motor control a and the top-down neural prediction g(μ), weighted by the neural inertial mass m z . Below, we appraise the BM prescribed by Eqs. (17) and (18) and note some significant aspects: (i) The derived RD suggests that both brain activities μ and their conjugate momenta p are dynamic variables. The instantaneous values of μ and p correspond to a point in the brain's perceptual phase space, and the continuous solution over a temporal horizon forms an optimal trajectory that minimizes the theoretical action, which represents sensory uncertainty. (ii) The canonical momentum p defined in Eq. (13) can be rewritten as p = m w (μ − f ) − ρ √ m w /m z ϕ . Accordingly, when the normalized correlation ρ is nonvanishing, the momentum quantifies combined errors in predicting changing states and sensory stimuli. Prediction errors propagate through the brain by obeying coupled dynamics according to Eqs. (17) and (18). (iii) Terms involving time-dependent ϕ in Eqs. (17) and (18) are identified as driving forces C i , i = μ, p, The sensory prediction error ϕ defined in Eq. (20) quantifies motor signals engaging the brain's nervous control in integrating the RD.
Equations (17) and (18) are the highlights of our formulation, which prescribe the brain's BM of semi-actively inferring the external causes of sensory inputs under the revised FEP. Note that the motor variable a is not explicitly included in our derived RD; instead, it implicitly induces nonautonomous sensory inputs ϕ(t) in the motor signal ϕ . The motor signal appears as a time-dependent driving force; accordingly, our Hamiltonian formulation bears a resemblance to the motor-control dynamics described by the Hamilton-Jacobi-Bellman (HJB) equation in the control theory (Todorov 2007). If one regards the Lagrangian Eq. (11) as a negative cost rate and the canonical momentum p as a costate, our IA is equivalent to the total cost function that generates the continuous-state HJB equations. In optimal control theory, the associated Hamiltonian function is further minimized with respect to the control signal, which we do not explicitly consider in this work. In our formulation, the motor signals are produced by the discrepancy between the sensory streams ϕ(t) and those predicted by the brain. The nonstationary data are presented to a sensorimotor receptor, whose field position in the environment is specified by the agent's locomotive motion. The neural observer continuously integrates the BM subject to a motor signal to perform the sensoryuncertainty minimization, thereby closing the perception and motion control within a reflex arc. When we neglect the correlation ρ between the sensory prediction modeled by Eq. (6) and the internal dynamics of predicting the neuronal state modeled by Eq. (5), we can recover the RD reported in the previous publication (Kim 2018), which demonstrates the consistency of our formulation.
In the present treatment, we consider only a single brain variable μ; accordingly, the ensuing BM specified by Eqs. (17) and (18) is described in a two-dimensional phase space. The extension of our formulation to the general case of the multivariate brain is possible by applying the same line of work proposed in (Kim (2018)). Under the independentparticle approximation, the multivariate Lagrangian takes the form where {μ} = (μ 1 , μ 2 , . . . , μ N ) denotes a row vector of N brain states that respond to multiple sensory inputs {ϕ} = (ϕ 1 , ϕ 2 , . . . , ϕ N ) in a general manner. Note that our proposed multivariate formulation is different from the state-space augmentation using the higher-order states [see Sect. 2.1]. In our case, multiple brain states are, for instance, the membrane potential, gating variables, and ionic concentrations, that can be viewed as the fluctuating variables on a coarsegrained time scale, influenced by Gaussian white noises [see Sect. 2.3]. Furthermore, the implication of our formulation in the hierarchical brain can be achieved in a straightforward manner as in (Kim 2018), which adopts the bidirectional facet in information flow of descending predictions and ascending prediction errors (Markov and Kennedy 2013; Michalareas et al. 2016). Note that in ensuing formulation, both descending predictions and ascending prediction errors will constitute the dynamical states governed by the closed-loop RD in the functional architecture of the brain's sensorimotor system. This feature is in contrast to the conventional implementation of the FEP, which delivers the backward prediction-belief propagation-as neural dynamics and the forward prediction error as an instant message passing without causal dynamics (Friston 2010a;Buckley et al. 2017).

Simple Bayesian-agent model: implicit motor control
In this section, we numerically demonstrate the utility of our formulation using an agent-based model, which is based on a previous publication (Buckley et al. 2017). Unlike that in the previous study, the current model does not employ generalized states and their motions; instead, the RD is specified using only position μ and its conjugate momentum p for incoming sensory data ϕ. Environmental objects invoking an agent's sensations can be either static or time dependent, and in turn, the time dependence can be either stationary (not moving on average) or nonstationary. According to the framework of active inference, the inference of static properties corresponds to passive perception without motor control a. Meanwhile, the inference of time-varying properties renders an agent's active perception of proprioceptive sensations by discharging motor signals ϕ via classic reflex arcs.
In the present simulation, the external hidden state ϑ is a point property, e.g., temperature or a salient visual feature, which varies with the field point x. As the simplest environmental map, we consider h(ϑ, a) = ϑ(x(a)) and assume that the sensory influx at the corresponding receptor is given by where z gp denotes the random fluctuation. The external property, e.g., temperature, is assumed to display a spatial profile as where ϑ 0 denotes the value at the field origin, and the desired environmental niche is situated at The biological agent that senses temperature is allowed to navigate through a one-dimensional environment by exploiting the hidden property. The agent initiates its motion from x(0), where the temperature does not accord with the desired value. In this case, the agent must fulfill its allostasis at the cost of biochemical energy by exploiting the environment based on where a(t) denotes a motor variable, e.g., agent's velocity. The nonstationary sensory data ϕ(t) are afferent at the receptor subject to noise z gp ; its time dependence is caused by the agent's own motion, i.e., ϕ(t) = ϑ(x(a(t))), which is assumed to be latent to the agent's brain in the current model. With the prescribed sensorimotor control, the rate of sensory data averaged over the noise is related to the control variable aṡ The neural observer is not aware of how sensory inputs at the proprioceptor are affected by the motor reflex control of the agent. In the case of saccadic motor control (Friston et al. 2012b), an agent may stand at a field point without changing its position; however, sampling the salient visual features of the environment through a fast eye movement a(t) makes the visual input nonstationary, i.e., ϕ(t) = ϑ(a(t)).
In Fig. 1, we depict streams of sensory data at the agent's receptor as a function of time. For this simulation, the latent motor variable in Eq. (25) is considered as which renders the agent's position in the environment as (0). For simplicity, we assume that this is hardwired in the agent's reflex pathway over evolutionary and developmental time scales. The figure shows that the agent, initially located at x(0) = 10, senses an undesirable stimulus ϑ(0) = 0.2; accordingly, it reacts by using motor control to determine an acceptable ambient niche. For this illustration, we assumed the environmental property at the origin to be ϑ 0 = 20. After a period of t = 5, the agent finds itself at the origin x = 0, where the environmental state is marked by the value ϑ = 20.
Having prescribed the nonstationary sensory data, we now set up the BM to be integrated by applying Eqs. (17) and (18) to the generative models below. We assume that the agent has already learned an optimal generative model; therefore, the agent retains prior expectations regarding the observations and dynamics. Here, for the demonstration, we consider the learned generative model in its simplest linear form Note that the motor control a is not included in the generative model, and the desired sensory data ϑ d , e.g., temperature, appear as the brain's prior belief of the hidden state. Accordingly, Eqs. (17) and (18) are reduced to a coupled set of differential equations for the brain variable μ and its conjugate momentum p aṡ Parameters α, β, and γ are proportional to the correlation ρ; see Eq. (19). Hence, they become zero when the neural response to the sensory inputs is uncorrelated with neural dynamics, which is not the case in general. Time-dependent driving terms appearing on the RHS of both equations, namely Eqs. (28) and (29), include the sensorimotor signal ϕ (μ; ϕ(t)) given in Eq. (20). The motor variable a, which drives the nonstationary inputs ϕ(t), is unknown to the neural observer in our implementation.
In the following, for a compact mathematical description, we denote the brain's perceptual state as a column vector Vector Ψ represents the brain's current expectation μ and the associated prediction error p with respect to the sensory causes, as encoded by the neuronal activities performed when encountering a sensory influx. Therefore, in terms of perceptual vector Ψ , Eqs. (28) and (28) are expressed as where relaxation matrix R is defined as and source vector S encompassing the sensory influx ϕ(t) is defined as Unless it is a pathological case, the steady-state (or equilibrium) solution ψ eq of Eq. (30) is uniquely obtained as We find it informative to consider the general solution Ψ (t) of Eq. (30) with respect to the fixed point ψ eq by setting To this end, we seek time-dependent solutions for the shifted measure ψ(t) as follows where δS = S(t) − S(∞). It is straightforward to integrate the above inhomogeneous differential equation to obtain a formal solution, which is given by Note that δS becomes zero identically for static sensory inputs; therefore, the relaxation admits simple homogeneous dynamics. In contrast, for time-varying sensory inputs, the inhomogeneous dynamics driven by the source term is expected to be predominant. However, on time scales longer than the sensory-influx saturation time τ , it can be shown that δS → 0; for instance, τ = 5 in Fig. 1. Therefore, for such a time scale, the inhomogeneous contribution in the relaxation diminishes even for time-varying sensory inputs, and the homogeneous contribution is dominant for further timedevelopment. The ensuing homogeneous relaxation can be expressed in terms of eigenvalues λ l and eigenvectors ξ (l) of the relaxation matrix R as where expansion coefficients c l are fixed by initial conditions ψ(0). The initial conditions ψ(0) represent a spontaneous or resting cognitive state. In Eq. (35), eigenvalues and eigenvectors are determined by the secular equation Then, the solution for the linear RD Eq. (30) is given by which is exact for perceptual inference, and legitimate for active inference on timescales t > τ. Before presenting the numerical outcome, we first inspect the nature of fixed points by analyzing the eigenvalues of the relaxation matrix R given in Eq. (31). First, it can be seen that the trace of R is zero, which indicates that the two eigenvalues have opposite signs, i.e., λ 1 = −λ 2 . Second, the determinant of R can be calculated as Therefore, if the correlation φ → 0, it can be conjectured that both eigenvalues are real. This is because Det(R) = λ 1 λ 2 → −1 − m z /m w < 0, which yields λ 2 1 = λ 2 2 > 0 using the first conjecture. Thus, we can conclude that the two eigenvalues are real and have opposite signs. Therefore, for φ = 0, the solution is unstable. In contrast, when the correlation is retained, Det(R) can be positive for a suitable choice of statistical parameters, namely m w , m z , and φ. In the latter case, the condition λ 1 λ 2 > 0 renders λ 2 l < 0 for both l = 1, 2. Accordingly, λ 1 and λ 2 that have opposite signs are purely imaginary, which makes the fixed point Ψ eq a center (Strogatz 2015). If we define λ 1,2 ≡ ±iω, the long-time solution of RD with respect to Ψ eq is expressed as where the state space is spanned by continuous brain state μ and its conjugate momentum p variables. The common fixed point is indicated by an orange bullet at the center of the orbits, which predicts the sensory cause incorrectly. [All curves are in arbitrary units.]

Fig. 3
Active perception: Time-development of the perceptual state inferring the external causes of sensory inputs altered by the agent's motor control. Blue and magenta curves depict the brain activity μ(t) and corresponding momentum p(t), respectively. In addition, the noisy curve indicates the nonstationary sensory inputs ϕ(t) entering the sensory receptor at instant t. For numerical illustration, we used σ w = 1.0, σ z = 10, and φ = −2.8. [All curves are in arbitrary units.] which specifies a limit cycle with angular frequency ω. Thus, according to our formulation, the effect of correlation on the brain's RD is not a subsidiary but a crucial component. Below, we consider numerical illustrations with finite correlation.
We exploited a wide range of parameters for numerically solving Eqs. (28) and (29) and found through numerical observation that there exists a narrow window in the statistical parameters σ w , σ z , and φ, within which a stable trajectory is allowed for a successful inference. This finding implies that the agent's brain must learn and hardwire this narrow parameter range over evolutionary and developmental timescales; namely, generative models are conditioned on an individual biological agent. We denote the instantaneous cognitive state as (μ(t), p(t)) for notational convenience.
In Fig. 2, we depict the numerical outcome from the perceptual inference of static sensory inputs. To obtain the results, we select a particular set of statistical parameters as σ w = 1.0, σ z = 10, and φ = −2.8, which specify the neural inertial masses m w = 4.6 and m z = 0.1 × m w and the coefficients that enter the RD, namely α = −0.60, β = −0.28, and γ = −2.8 (κ = 10).
In Fig. 2Left, we depict the brain variable μ as a function of time, which represents the cognitive expectation of a registered sensory input under the generative model [Eq. (26)] for three values, namely ϕ = 10, 15, and 20. For all illustrations, the agent's prior belief with regard to the sensory input is set as which is indicated by the horizontal dotted line. The blue curve represents the case in which sensory data are in line with the belief. The RD of the perceptual inference delivers an exact output (μ eq , p eq ) = (20, 0), where μ eq and p eq are the perceptual outcome of the sensory cause and its prediction error, respectively. Note that μ eq and p eq correspond to the temporal averages of μ(t) and p(t), respectively. The other two inferences underscore the correct answer. Figure 2Right corresponds to the case of a single sensory data ϕ = 4.0, which the standing agent senses at the field point x = 2. The ensuing trajectories from all three initial spontaneous states have their limit cycles in the state space defined by μ and p. We numerically determined the fixed point to be (μ eq , p eq ) = (−65.6, −306) and the two eigenvalues of the relaxation matrix R to be (λ 1 , λ 2 ) = (1.84i, −1.84i), which are purely imaginary and have opposite signs. Again, the perceptual outcome does not accord with the sensory input; it deviates significantly.
Next, in Fig. 3, we depict the results for active inference, which were calculated using the same generative parameters used in Fig. 2. The agent is initially situated at x(0) = 2, where it senses the sensory influx ϑ(0) = 4, which does not match the desired value ϑ d = 20. Therefore, the agent reacts to identify a comfortable environmental niche matching its prior belief, which generates nonstationary sensory inputs at the receptors (Fig. 1). The brain variable μ initially undergoes a transient period at t ≤ 5. The RD commences from the resting condition (μ(0), p(0)) = (0, 0) and then develops a stationary evolution. Furthermore, we numerically confirmed that the brain's stationary prediction μ eq , which is the brain's perceptual outcome of the sensory cause, is close to but not in line with the prior belief ϑ d . The stationary value p eq is estimated to be approximately 8.0, which is the average of the stationary oscillation of prediction error p(t).
In Fig. 4, the trajectory corresponding to that in Fig. 3 is illustrated in blue in the perceptual state space spanned by μ and p, including two other time developments from different choices of initial conditions. All data were calculated using the same generative parameters and sensory inputs used for Fig. 3. Regardless of the initial conditions, after each transient period, the trajectories approach stationary limit cycles about a common fixed point, as seen in the case of static sensory inputs in Fig. 2. The fixed point Ψ eq and stationary frequency ω of the limit cycles are not affected by initial conditions, which are solely determined by the generative parameters m w , m z , and φ and the prior belief ϑ d for a given sensory input ϕ [Eqs. (33) and (36)]. In addition, we numerically observed that the precise location of the fixed points is stochastic, thereby reflecting the noise from the nonstationary sensory influx ϕ.
In the framework of active inference, motor behavior is attributed to the inference of the causes of proprioceptive sensations (Adams et al. 2013), and in turn, the prediction errors convey the motor signals in the closed-loop dynamics of perception and motor control. In Fig. 5, we depict the sensorimotor signals ϕ (μ; ϕ(t)) that appear as timedependent driving terms in Eqs. (28) and (29). In both figures, the agent is assumed to be initially situated such that it can sense the sensory data ϕ(0) = 4. After an initial transient period elapses, the motor signals exhibit a stationary oscillation about average zero in Fig. 5Left, implying the successful fulfillment of the active inference of nonstationary sensory influx matching the desired belief ϑ d = 20. The amplitude of the motor signal shown by the blue curve is smaller than that shown by the red curve, which is also reflected in the size of the corresponding limit cycles in Fig. 4. The prediction-error signal from the plain perception exhibits an oscillatory feature in the gray curve, which arises from the stationary time dependence of the brain variable μ(t). The amplitude shows a large variation caused by the significant discrepancy between the static sensory input ϕ = 4 and its prior belief ϑ d = 20. In Fig. 5Right, we repeated the calculation with another value: ϑ d = 10. In this case, the prior belief ϑ d regarding the sensory input does not accord with stationary sensory streams. Therefore, the blue and red signals for active inference oscillate about the negatively shifted values from average zero. In contrast to Fig. 5Left, the error-signal amplitude of the static input is reduced because the difference between the sensory data and prior belief decreases.
Next, we consider the role of correlation φ in the brain's RD, whose value is limited by the constraint |φ| ≤ √ σ w σ z . To this end, we select three values of φ for the fixed variances σ w and σ z , and we integrate the RD for active inference. In Fig. 6, we present the resulting time evolution of the brain states μ for the initial condition (μ(0), p(0)) = (0, 0). In this figure, the conjugate momentum variables are not shown. The noticeable features in the results include the changes in Motor signals ϕ (μ; ϕ(t)) as a function of time t evoked by the discrepancy between the nonstationary sensory stream and its top-down prediction [Eq. (20)]. Here, we set the prior belief ϑ d = 20 (Left) and ϑ d = 10 (Right). The blue and red curves represent the results from the initial condition (μ(0), p(0)) = (0, 0) and (0, 100), respectively.
The gray curves represent the corresponding signals from the plain perception of the static sensory input. All data were obtained by setting the statistical parameters as σ w = 1.0, σ z = 10, and φ = −2.8. [All curves are in arbitrary units.] the fixed point and the amplitude of the stationary oscillation with correlation. The average value of μ(t) in the periodic oscillation corresponds to the perceptual outcome μ eq of the sensory data in the stationary limit. We remark that for all numerical data presented in this work, we selected only negative values for φ. This choice was made because our numerical inspection revealed that positive correlation does not yield stable solutions.
In Fig. 7, as the final numerical manifestation, we show the temporal buildup of the limit cycles in the perceptual phase space; however, this time, we fix σ w while varying σ z and φ. To generate the red, blue, and gray curves, the tuning parameter κ was selected as κ = 50, 10, and 100, respectively. The resulting fixed points are located approximately at the center of each limit cycle, which are not shown. Similar to that in Fig. 6, it can be observed that the positions of the fixed point and amplitudes of oscillation are altered by variations in the statistical parameters. Evidently, a different set of parameters, namely σ w , σ z , and φ, which are the learning parameters encoded by the brain, result in a distinctive BM of active inference.
Here, we summarize the major findings from the application of our formulation to a simple nonstationary model. The brain's BM, i.e., Eqs. (28) and (29), employs linear generative models given in Eqs. (26) and (27).
(i) The steady-state solutions of the RD turn out to be a center about which stationary limit cycles (periodic oscillations) are formed as an attractor (Friston and Ao 2012a) in the perceptual phase space, which constitute the brain's nonequilibrium resting states. (ii) The nonequilibrium stationarity stems from the pair of purely imaginary eigenvalues of the relaxation matrix with opposite signs, given by Eq. (31); the equal magnitude specifies the angular frequency ω of the periodic trajectory. Fig. 6 Time evolution of the brain variable μ. Here, we vary the correlation φ for fixed variances σ w = 1.0 and σ z = 10. The red, blue, and gray curves correspond to φ = −3.0, −2.8, and −2.6, respectively. For all data, the agent is environmentally situated at x = 2, where it senses the transient sensory inputs ϕ(t) induced by the motor reflexes at the proprioceptive level. The agent's initial cognitive state is assumed to be (μ(0), p(0)) = (0, 0), and the prior belief is set as ϑ d = 20. [All curves are in arbitrary units.] (iii) Centers are determined by generative parameters and the prior belief for a given sensory input [Eq. (33)], which represents the outcome of active inference and the entailed prediction error. (iv) The theoretical assumption of the statistical dependence of two generative noises describing the brain's expectation of the external dynamics and sensory generation is consequential to ensuring a stationary solution. Furthermore, based on numerical experience, a negative covariance is necessary for obtaining stable solutions using the current model.

Fig. 7
Limit cycles in the perceptual phase space spanned by the brain state μ and its conjugate momentum p. Here, we considered several sets of σ z and φ for a fixed σ w = 1.0. The red, blue, and gray curves were obtained from (σ z , φ) = (50, −6.6), (10, −2.8), and (100, −9.5), respectively. For all data, the agent's initial cognitive state is assumed to be (μ(0), p(0)) = (0, 0), and the prior belief is set as ϑ d = 20. The agent is environmentally situated at x = 2, where it senses the transient sensory inputs ϕ(t) induced by the motor reflexes at the proprioceptive level. [All curves are in arbitrary units.]

Concluding remarks
In the present study, we continued our effort to make the FEP a more physically principled formalism based on our previous publication (Kim 2018). We implemented the FEP in the scope of the principle of least action by casting the minimization scheme to the BM described by the effective Hamiltonian equations in the neural phase space. We deconstructed some of the theoretical details in the first part, which are embedded in the formulation of the FEP, while comparing our approach with other currently prevailing approaches. In the second half, we demonstrated our proposed continuousstate RD in the Bayesian brain using a simple model, which is biologically relevant to sensorimotor regulation such as motor reflex arcs or saccadic eye movement. In our theory, the time integral of the induced IFE in the brain, not the instant variational IFE, performs as an objective function. In other words, our minimization scheme searches for the tight bound on the sensory uncertainty (average surprisal) and not the instant sensory surprisal.
To present the novel aspects of our formulation, this study focused on the perceptual inference of nonstationary sensory influx at the interface. The nonstationary sensory inputs were assumed to be unknown or contingent to the neural observer without explicitly engaging in motor-inference dynamics in the BM. Instead, we considered that the motor signals are triggered by the discrepancies between the sensory inputs at the proprioceptive level and their top-down predictions. They appeared as nonautonomous source terms in the derived BM, thus completing the sensorimotor dynamics via reflex arcs or oculomotor dynamics of sampling visual stimuli. This closed-loop dynamics contrasts with the gradient-descent implementation, which involves the double optimization of the top-down belief propagation and the motor inference in message-passing algorithms. In our present formulation, the sensorimotor inference was not included; however, a mechanism of motor inference can be included explicitly by considering a Langevin equation for a sensorimotor state. This procedure extends the probabilistic generative model by accommodating the prior density for motor planning for active perception, which is similar to what was done in Bogacz (2020).
By integrating the Bayesian equations of motion for the considered parsimonious model, we manifested transient limit cycles in the neural phase space, which numerically illustrate the brain's perceptual trajectories performing active perception of the causes of nonstationary sensory stimuli. Moreover, we revealed that ensuing trajectories and fixed points are affected by the input values of the learning parameters (both diagonal and off-diagonal elements of the covariance matrix) and prior belief regarding sensory data. The idea of exploring the effect of noise covariance was purely from the theoretical insight without a supporting empirical evidence, which allowed us to drive a stable solution in perceptual and motor-control dynamics. We did not attempt to explicate in detail the effect of neural inertial masses (precisions) and correlation (noise covariance) on the numerically observed limit cycles. This was because of the numerical limitation set by the presented model, which permits stable solutions in a significantly narrow window of statistical parameters. In neurosciences, it is commonly recognized that neural system dynamics implement cognitive processes influencing psychiatric states (Durstewitz et al. 2020). We hope that the key features of our manifestation will serve to motivate and guide further investigations on more realistic generative models with neurobiological and psychological implications.
Finally, we mention the recent research efforts on synthesizing perception, motor control, and decision making within the FEP (Friston et al. 2015Biehl et al. 2018;Parr and Friston 2019;van de Laar and de Vries 2019;Tschantz et al. 2020;Da Costa et al. 2020a). The underlying idea of these studies is rooted in machine learning (Sutton and Barto 1998) and the intuition from nonequilibrium thermodynamics (Parr et al. 2020;Friston 2019), and they attempt to widen the scope of active inference by incorporating prior beliefs regarding behavioral policies. The new trend supplements the instant IFE to the future expected IFE in a time series, and it formulates the adaptive decision-making processes in action-oriented models. The assimilation of this feature needs to be studied in depth (Millidge et al. 2020b;Tschantz et al. 2020). We are currently considering a formulation of motor inference together with the assimilation of extended IFEs in the scope of the least action principle.
Funding Not applicable.

Conflict of interest Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.