$\mathrm{\Phi }$ was calculated with autoregressive $\mathrm{\Phi }$ algorithm (${\mathrm{\Phi }}_{\text{AR}}$) [28]. This approach was developed for performing calculations on time-series data, because the estimation of entropies in the definition (1) requires to estimate the whole distribution of data, to gather statistics about $2^{\text{N}}$ states of a system of $N$neurons. For real time-series data, this distribution is often underestimated and yields to instability of $\mathrm{\Phi }$ [28]. To build binary time-series data for calculations, a timeline was divided into equal bins. If there was a spike in the bin, this bin has value 1, and in other case it has value 0. A bin size is used as a compromise, in order to keep all patterns of neuron spiking and avoid excessive information. Such binning method was used for calculation of entropies, for example in [29].

This algorithm has a parameter $\mathrm{\Delta }t$, which means how much of sequential states are considered (for how many steps the information is generated). This is a temporal scale for a system, it shows a degree to which an information in a network predicts a future network state given an earlier state separated in time by a lag $\mathrm{\Delta }t$ [30]. So, $\mathrm{\Delta }t$ was determined experimentally by calculating $\mathrm{\Phi }$ with different values of $\mathrm{\Delta }t$ and finding such $\mathrm{\Delta }t$ which gives a maximal average value of $\mathrm{\Phi }$ for all learning days.

This classical definition (1) cannot be used for calculation of $\mathrm{\Phi }$. It is based on estimation of a partition minimizing $\mathrm{\Phi }$ , which can only be found by a brute-force search. A number of possible ways to partition a system of N elements is the Bell number $B_N$ [31]. As a system size N increases, this number grows faster than a factorial. Thus, another calculation method, so-called autoregressive $\mathrm{\Phi }$ (${\mathrm{\Phi }}_{\text{AR}}$) was used [28]:

(2)$${\mathrm{\Phi }}_{\mathrm{AR}}=\frac{\underset{M}{\min }\left(\frac{1}{2}\ln \frac{\det (\mathrm{\Sigma }(X))}{\det \left(\mathrm{\Sigma }\left(E^X\right)\right)}-{\sum }_{k=1}^2\frac{1}{2}\ln \frac{\det \left(\mathrm{\Sigma }\left(M_k\right)\right)}{\det \left(\mathrm{\Sigma }\left(E^{M_k}\right)\right)}\right)}{\frac{1}{2}\ln \left[\underset{k}{\min }\left\{{(2\pi e)}^{\left|M_k\right|}\det \left(\mathrm{\Sigma }\left(M_k\right)\right)\right\}\right]}$$

where $\mathrm{\Sigma }\left(X\right)$ is the covariation matrix of variable $X;$$E^X$ is the autoregression residuals (errors of regression that predicts value of $X$ at time moment $t+\mathrm{\Delta }t$ based on value of $X$ at time moment $t;$$M_{k,}k\in \{0,\mathrm{\hspace{0.25em}1}\}$ are the two system subsets, and $\left|M_k\right|$ is a size of a subset $M_k$. Subsets are selected in a way to minimize the resulting ${\mathrm{\Phi }}_{\text{AR}}$ value. Only two subsets are selected in this method, the number of candidate partitions is defined by a Stirling number of the second kind and grows exponentially. But for a large number of neurons a calculation in a considerable time is still impossible. Thus, an approximation was applied. Experimentally we found that the limit of 15 neurons is the point when computations become intractable. To perform calculations, at each session we selected 15 neurons with the most variance of activity, as it was done in [26]. In previous works, other approximations based on simplified partition were used [32].

In some cases, this algorithm could not return a value. It is known [26], that it happens when some neurons had little or no variance in their activities — bins were almost always on or almost always off. In this case, the one neuron with the least variance was excluded from analysis and calculation was performed again, and so on. If an algorithm was unable to calculate ${\mathrm{\Phi }}_{\text{AR}}$ on any neurons in some period, this period was excluded from statistical analysis. Each obtained ${\mathrm{\Phi }}_{\text{AR}}$ value was normalized by a number of neurons used in this calculation to make it possible to compare obtained values and investigate its dynamics through each animal learning sessions.

To achieve a higher time resolution, each session was divided into 8 equal periods and information integration coefficient ${\mathrm{\Phi }}_{\text{AR}}$ was also calculated for each period. Calculations were performed in MATLAB. Neural and behavioral data manipulations were performed with FMAToolbox open-source library (http://fmatoolbox.sourceforge.net) [33].

We used an open dataset [34, 35] in order to analyze behavior and neuronal activity during learning, so details of the experimental design may be found in [35]. Animal experiments were approved by the Institutional Animal Care and Use Committee (IACUC) at New York University Medical Center, as it was mentioned in an original manuscript [35]. In short, male Long-Evans rats (N = 3, 300 g, 3 months) were initially trained a simple spatial task in a linear track. During behavior sessions the animals ran from one end of the track to the other to get water from a water well. After shaping or initial training sessions, microelectrode arrays were implanted into the amygdala and the dorsal hippocampus. After a rehabilitation period of 5 days, animals started learning the second procedure. Two air-puff sites were located at the equal distance from the middle of the track and applied at one of two sites when the rat moved in a particular direction. The direction and the site varied pseudorandomly across training days.

Along with information integration coefficient we quantified learning success for every period of each session. Learning success was assessed as a number of rewards in a period divided by a length of a period. Correlations between the two variables in each period were calculated using Spearman’s correlation; we also calculated Spearman’s $\rho$ using Student’s t test. The Bonferroni correction for multiple comparisons was performed.