The mirror system Let us begin with a simple dissipative dynamical system tracking the evolution of a dependent variable \(x\left(t\right)\) that decays exponentially in time:

where \(a\) is a positive constant.

In addition, we describe the environment in terms of a mirror system that evolves according to a new variable \(\widetilde\left(t\right)\), where this system gains energy at the same rate as the original system in Eq. (1) loses energy, such that:

$$\dot}=a\widetilde$$

(2)

The need for this mirror system approach can be motivated by the limitations of alternative approaches, as outlined in Appendix 1.

We can then cast Eqs. (1) and (2) in the form of the following Lagrangian \(\mathcal\):

$$\mathcal=\frac\left(\dot\widetilde-x\dot}\right)+ax\widetilde$$

(3)

which we can verify by recovering the equations of motion via the Euler–Lagrange equations:

$$\begin \frac}}}}} - \frac}\left( }}}}}}} \right) = 0~~~~ \Rightarrow ~~~~\dot = - ax \hfill \\ \frac}}}} - \frac}\left( }}}}}} \right) = 0~~ \Rightarrow ~~~\dot} = a\tilde \hfill \\ \end$$

(4)

The energyFootnote 1\(E\) associated with the Lagrangian in Eq. (3) is then obtained via the Legendre transformation:

$$\begin E & = \frac}}}}}}\dot} + \frac}}}}}\dot - } \\ & = - ax\tilde \\ \end$$

(5)

The solutions of Eqs. (1) and (2) read:

$$\begin x\left( t \right) & = c_ e^ \\ \tilde\left( t \right) & = c_ e^ \\ \end$$

(6)

where \(_\) and \(_\) are integration constants.

Using Eqs. (5) and (6) we then find that the energy is conserved, as required (Fig. 2)

$$E = - ac_ c_ \Rightarrow \dot = 0$$

(7)

For chosen values of \(a = - 0.4\), \(c_ = 2\), \(c_ = 4\), we show time (arbitrary units) on the horizontal axis and the dependent (dep.) variables on the vertical axis—consisting of the time course of \(x\left( t \right)\) (solid black line), the corresponding mirror time course \(\tilde\left( t \right)\) (dashed line) from Eq. (6), and the resultant constant energy \(E\) (red line) time course from Eq. (7)

Bayesian model inversion We use dynamic causal modelling (DCM) with the statistical parametric mapping (SPM) software to extract estimates of model parameters (Friston et al. 2003) for a given time course. This routine sets all free parameters (such as coupling strengths) to Bayesian priors of zero—providing starting points from which the model inversion then searches for Bayesian posterior values that best explain the underlying data, by using variational Laplace to estimate the variance of states (Roebroeck et al. 2011).

Specifically, we begin by setting the model states in terms of the initial conditions of the dependent variable at \(t=0\). We also set the prior means of model parameters to zero and set prior variances of unity, thereby allowing the posteriors of the model parameters to deviate from the prior means. We set up an observation function consisting of the dependent variable and apply an external driving input that is initialized as a randomized time course with mean zero. Furthermore, we set the precisions of observation noise, state noise, and of exogenous causes and run Bayesian model inversion (spm_LAP in spm12) to infer the best fits of the model parameters and external driving input for a given dataset. This routine returns a summary statistic known as the variational free energy, which represents a trade-off between the accuracy and complexity of a given model. In other words, for a given level of accuracy, model evidence is penalized for each additional degree of freedom used, thereby avoiding an potential overfitting problem. Full details of this routine are provided in the accompanying code.

Resting-state fMRI We next use empirical neuroimaging timeseries to obtain estimates of energy across a 100-region brain parcellation (Schaefer atlas) of the cortex (Schaefer et al. 2018).

Regional fMRI blood-oxygen-level-dependent (BOLD) resting-state data were sourced from the 1200-subject release of the Human Connectome Project (HCP) (Essen et al. 2013). The acquisition of these data was performed using a Siemens Skyra 3 Tesla MRI scanner with a 3T imaging capability. The scanning parameters included a repetition time (TR) of 720 ms, an echo time (TE) of 33 ms, a flip angle of 52 degrees, and a voxel size of 2 mm isotropic. The scan covered 72 slices with a matrix size of 104 × 90 and a field of view (FOV) of 208 × 180 mm, employing a multiband acceleration factor of 8. From this dataset, 100 subjects aged between 22 and 35 were selected at random, with three individuals excluded due to missing data.

To mitigate the effects of global movement and respiratory-induced artifacts, we ran the fMRI data through noise reduction via the ICA-FIX based pipeline. FMRIB's ICA-based Xnoiseifier (ICA-FIX) serves as a preprocessing framework aimed at cleaning fMRI data by eliminating components associated with noise. This framework integrates two distinct techniques: independent component analysis (ICA) and FMRIB's Automated Removal of Motion Artefacts (FMRIB-AROMA) (Griffanti et al. 2014). We then used the 100-region Schaefer cortical atlas in MNI152 2 mm standard space to extract the mean regional fMRI signal for each subject.

Neuromaps Finally, we assessed whether our empirical resting-state fMRI-derived regional LTI energy estimates corresponded to any pre-existing regional measures of biologically-based notions of energy. To this end, we used the Neuromaps toolbox to obtain structural, functional, metabolic and electrophysiological maps, parcellated according to the same 100-region Schaefer cortical atlas as used with the resting-state fMRI datasets. These maps include regional cerebral metabolic rate (rCMRGlu), derived from 18-flurodeoxyglucose positron emission tomography (FDG-PET)—frequently used to estimate metabolic energy consumption in vivo. The maps also include regional cerebral blood flow (CBF) derived from MRI arterial spin labelling (ASL), which is tightly coupled to rCMRGlu in the healthy resting brain.

The list of 16 Neuromaps we used is as follows: cerebral blood flow, blood volume, metabolic rate, and metabolic rate of glucose (Vaishnavi et al. 2010); the first three diffusion map embedding gradients of group-averaged functional connectivities (Margulies et al. 2016); MEG timescales, together with alpha, beta, delta, low gamma, high gamma, and theta power distributions (Essen et al. 2013; Shafiei et al. 2022); MRI cortical thicknesses and T1w/T2w ratios (Glasser et al. 2016).

Significance testing To account for the contribution of spatial autocorrelation between maps, we correlated our subject-averaged LTI energy map with \(N=\text\) spatial null maps for significance testing (Markello and Misic 2021). This allows us to calculate a p-value by assessing where a given correlation strength lies within a sorted list of \(\text\) correlation strengths, each of which is calculated from maps that have a different spatially randomized (spinpermuted) segmentation. One such significance test was performed for each of the 16 Neuromaps.

As a further test of significance, we repeated all correlations with mean signal values, rather than with mean energy values. This is done in order to assess whether a given result, once spatial autocorrelation is accounted for, is specifically due to the form of the LTI energy expression (see Eq. (15) in Results), as opposed to some more basic property of the time course that can be captured just by its mean value.