Causal relationship between local field potential and intrinsic optical signal in epileptiform activity in vitro

The directed causal relationship were examined between the local field potential (LFP) and the intrinsic optical signal (IOS) during induced epileptiform activity in in vitro cortical slices by the convergent cross-mapping causality analysis method. Two components of the IOS signal have been distinguished: a faster, activity dependent component (IOSh) which changes its sign between transmitted and reflected measurement, thus it is related to the reflectance or the scattering of the tissue and a slower component (IOSl), which is negative in both cases, thus it is resulted by the increase of the absorption of the tissue. We have found a strong, unidirectional, delayed causal effect from LFP to IOSh with 0.5-1s delay, without signs of feedback from the IOSh to the LFP, while the correlation was small and the peaks of the cross correlation function did not reflect the actual causal dependency. Based on these observations, a model has been set up to describe the dependency of the IOSh on the LFP power and IOSh was reconstructed, based on the LFP signal. This study demonstrates that causality analysis can lead to better understanding the physiological interactions, even in case of two data series with drastically different time scales.


Results
Three methods have been applied to evoke the epileptiform activity in rat neocortical slice preparations: electric stimulation (STIM, 5 slices), application of Mg 2+ -free artificial cerebrospinal fluid (4 slices) and 4-aminopiridine (4AP) treatment (4 slices). All the three methods evoked epiletiform activity by means of recurrent population burst discharges in the slices. There were specific differences between the burst discharges which were characteristic for the method of induction and there were large variation among slices as well. The induction-specific characteristic properties include the amplitude, the intraburst and interburst frequencies as well as the temporal development of the activity during the 1 hour long measurements 14 . The IOS was measured parallel, by quantifying the relative changes in the reflected or transmitted light flux over the somatosensory cortex. In 12 slices, the IOS signal was measured in interface setup, which allowed the measurement of the reflected light, but in one slice, the Mg 2+ -free induced activity was measured in submerged setup, which made possible the measurement of the transmitted light. An example of the Mg 2+ -free evoked epileptiform activity, the parallel recorded LFP and reflected IOS are shown in Fig. 1B. The slow and the fast components of the IOS signal. The general temporal development of the reflected IOS signal during the epileptic activity induction is consists of two components: The first component Time delayed causality between LFP power and the fast component of IOS. The sampling frequency of the LFP signal was 1 kHz, while the sampling rate of the parallel IOS recordings were only 2 Hz. As the dynamics of IOS is much slower than the LFP, this low sampling rate was enough to track the temporal changes of the IOS. However, the causality analysis requires fully synchronous time series. To reach this, the LFP signal was undersampled, by calculating the sum of the squared LFP amplitudes (the LFP power) within each 0.5 s long sampling interval.
The CCM method of Sugihara 12 , extended for delayed effects by Ye et al. 13 , was implemented in Scilab and applied to reveal the causality between the 1 hour long LFP and the IOSh recordings. Applying the time-delayed CCM analysis to the downsampled LFP power and to the IOSh recordings from each slice in the reflected light setup, we constructed the cross-map functions, which express the strength of the causal connections in both directions, as a function of the time delay between the two time series (Fig. 3). The resulting cross-map functions were averaged for slices sharing the same induction process: electric stimulation, Mg 2+ -free solution and 4AP treatment.
The causality analysis revealed similar structure in all the three cases. The cross-map functions have two significant peaks: a sharper peak for the LFP → IOSh causal direction, where the LFP precedes the IOSh with 1 s in case of electric stimulation and with 0.5 s in Mg 2+ -free and 4AP cases (Fig. 3A,C,E black lines, peaks on the www.nature.com/scientificreports www.nature.com/scientificreports/ negative half of the x axis) and a wider peak in the IOS → LFP direction, where IOSh follows LFP with 2.5 s time delay in the electric stimulation case and 1.5 s delay in the Mg 2+ -free and 4AP cases (Fig. 3A,C,E red lines, peaks on the positive half of the x axis).
In all cases, the peaks of the IOSh → LFP causality function are located at an anti-causal time delay, which means, that the presumed cause, the IOSh, in these cases, follows the presumed effect (the LFP) in time. Assuming that there is no difference in the observational delay, this is clearly impossible if the causal premise is correct. But, how can these anti-causal peaks be interpreted? Sugihara's CCM search for the fingerprints of the cause in the caused time series, thus a delayed CCM peak, means that the past cause can be reconstructed best from a later section of the caused time series. In our case, the real cause, the past LFP, can be reconstructed from the IOSh 0.5-1 s later. In the reverse (IOSh → LFP) direction, the algorithm checks, if the IOSh can be determined from the LFP. Thus the anti-causal peaks in the IOSh → LFP direction mean that the future of the IOSh can be determined based on the earlier LFP. We expect the appearance of such an anti-causal peak if the caused dynamics are more deterministic and non-chaotic, thus their future can be determined from the cause for a long time.
Thus, the conclusion of the causality analysis is that the LFP drives the IOSh activity, which follows quite predictably the LFP dynamics. There were no significant peaks on the negative half axis of the cross-map functions in the IOSh → LFP direction, which means that we did not find evidence for a feedback effect from the IOSh to the LFP in the examined range of time delays.
As the cross-map functions showed similar structures in all the three cases, we can assume that although the mechanisms generating the epileptic activity are different in the three cases, the causal link, thus the underlying mechanism, connecting the LFP to the IOSh may be similar.
Theoretically, the cross-map between two fully deterministic systems should reach the amplitude 1 for data series long enough. Thus, maxima smaller than 1 could be the effect of noise (which is non-deterministic) and www.nature.com/scientificreports www.nature.com/scientificreports/ that the recordings are not long enough. We found, that the amplitude of the causality peak was highest for the electrical stimulation case (mean 0.62 ± 0.09 SE), smaller for the Mg 2+ -free (mean 0.44 ± 0.16 SE) and smallest for the 4AP evoked activity (mean 0.3 ± 0.11 SE). The difference within the amplitudes could be the result of the different activity level (less discharges) within the cases as well as the different signal to noise ratio.

Causality versus correlation
We have compared the cross-map function to the traditional cross-correlation function. The cross-correlation function between LFP power and IOSh also showed similar structure in all cases, albeit with some differences (Fig. 3B,D,F blue lines). All the three correlation functions had a significant negative peak on the negative half axis (corresponding to the LFP lead) although at different delays: at −4.5 s in case of electric stimulation (mean −0.36 ± 0.007 SE), at −2.5 s in the Mg 2+ -free case (mean −0.21 ± 0.06 SE) and at −1 s in the 4AP case (mean −0.16 ± 0.016 SE). The cross-correlation function showed a second significant but positive peak on the positive half axis only at the Mg 2+ -free case (mean −0.11 ± 0.004 SE), but increased monotonically in the two other cases. Considering, that for very large time delays, the two data series become independent, thus the cross-correlation function should decay to zero, this monotonic increase implies the existence of a second even wider positive peak shifted far into the positive direction in these cases as well.
Comparison of the cross-correlation functions (Fig. 3B,D,F blue lines) to the LFP → IOSh causality function ( Fig. 3A,C,E black lines) shows, that none of peaks of the cross-correlation functions corresponds to the real delay of the causal effect represented by the peaks of the cross-map functions. We concluded that the real causal effect is not clearly reflected in the peaks of the cross-correlation functions in these cases.
As the short delay (δ < 0.5 s −1 s) causal effect from the LFP to the IOSh is not accompanied by high correlation with similar short delay, the effect from the LFP to the IOSh should take a form, which can not be well approximated with the linear dependence. One possibility is that, instead of the IOSh, the derivative of the IOSh depends on the LFP power. This assumption is supported by the shape of the cross-correlogram: highest steepness is close to zero time delay.
Thus, we also calculated the cross-correlation function between the LFP power and the time derivative of the IOSh signal. We found, that the peaks of the causality (Fig. 3A,C,E black lines) well matches to the negative peaks of the LFP-dIOSh/dt cross-correlation function (Fig. 3B,D,F black lines) in all the three cases. This implies the possibility, that the temporal derivative of the IOSh depends on the LFP power during the epileptiform discharges in the tissue.

Autonomous dynamics of IOSh without discharges
According to the above results, a causal link clearly exists between LFP and IOSh. However, between the burst discharges, the IOSh decays towards the baseline, which, in lack of observable electric activity, seems to be a result of autonomous dynamics (Figs 1E and 2). In order to determine the form of the decay functions, we collected the IOSh time series during those interburst periods which were longer than 25 s from all slices. The amplitude of the IOSh was normalized to 1 for each period and averaged over all bursts within a slice. The logarithm of the absolute value of the IOSh showed clear linear temporal dependence, during large part of these interburst periods in all cases, which is the hallmark of exponential decay dynamics. Thus, a linear function was fitted to the logarithm of the mean IOSh values between 2.5 and 25 s for each slices. Finally the mean and the SD of the decay time constant were calculated for the three forms of the epileptic activity. The mean decay time constant was the shortest in the Mg 2+ -free evoked activity with τ 1 = 30.3 ± 3.2 s (11, 9 and 11 intervals from 3 slices); the decay was longer in the 4AP treated slices: τ 1 = 54.9 ± 18.7 s (15, 7, 7, 11 intervals from 4 slices) and the longest for the electrically stimulated slices: τ 1 = 73 ± 19.  www.nature.com/scientificreports www.nature.com/scientificreports/

Reconstruction of IOSh based on LFP
By combining these observations, a simple formula is inferred, approximately describing the dependency of the IOSh on the LFP power in the form of a differential equation: where W(t) scales the causal effect from LFP to IOSh. Here we assumed, that the decay constant τ 1 is constant through whole recordings. This simple model has been used to estimating the IOSh signal based on the known LFP power in the Mg 2+ -free and 4AP elicited cases, where the epileptiform activity has developed autonomously, without external stimulation. The reconstructed IOSh were filtered similarly by subtracting the moving window average as the real observed signal and were compared to it. The optimal model parameters were determined by numerical optimization, minimizing the Root Mean Square Error (RMSE) of the IOSh reconstruction. First, W(t) was set to a constant value according to the negative peak of the correlation between LFP power and the time derivative of the IOSh, but we found, that the residual error of the IOSh reconstruction showed clear temporal tendency in majority of the recordings. Thus we introduced a time dependent weight factor W(t) in form of an exponential decay: The optimal W 0 and the τ 2 time constant was determined by a grid search. The optimization converged to a minimum of the RMSE in 6 cases out of the 7 slices (2 out of 3 Mg 2+ -free and all the 4 4AP slices). Thus we concluded, that the exponentially decaying effective connection strength was a reasonable description of the long term temporal development in the majority of the cases. The optimal time constant were found to be τ 2 = 2839 ± 582 s (mean and SD) for the Mg 2+ -free cases and τ 2 = 2061 ± 1121 s (mean and SD) for the 4AP cases.
The comparisons between the reconstructed and the observed IOSh show that reconstruction follows closely the actual observations during the 1 h long experiment, through several dynamical changes in both the Mg 2+ -free and the 4AP cases. An example for Mg 2+ -free elicited activity is shown in Fig. 5 and another for 4AP in Fig. 6. Note, the large difference between the two dynamics and waveshape.
Simulations. Based on the derived formula for the LFP → IOSh dependence, simulations have been run in order to cross-validate the causality analysis results. www.nature.com/scientificreports www.nature.com/scientificreports/ The simulations consist of a driving variable (X), a logistic chaotic oscillator, from which a spike-like discharge activity is derived, by raising it to the 4 th power, and a second driven variable exhibiting linear dynamics and exponential decay corresponding to the inferred model of the IOSh: The time delay δ of the driver effect have been varied and the causal relationship have been measured by the cross-map function. These simulations showed, that the exact value of the time delay between the two variables was inferred precisely by the peaks of the cross-map function, and even the shape of the cross-map function and the anti causal peak in the IOSh → LFP direction resembles very much to the observed one (Fig. 7). Based on these simulations, we can conclude that the time delay of the causal effect can be determined precisely by this method, in this case.

Discussion
Two components of the IOS signal have been distinguished during induced epileptic activity in in vitro cortical slices. They were different not only in their time scale, but presumably the underlying mechanisms as well. The faster, activity-dependent component (IOSh) was positive in transmitted light and negative in reflected light measurements. It can be interpreted as the decrease of the scattering of the tissue, caused by swelling of the cells due to the activation and underlying movement of ions across the membranes followed by water 15 . In this processes not only the neurons but also the glial cell might play significant role 16,17 . The later types of cells are more prone to swell, but the glial reaction is preceded by neuronal activation 18 . During this activation, an abundance of excitatory transmitters are released and intensive ionic movements take place increasing significantly the extracellular K + concentration, which is then buffered in glial cell and causes its swelling 19 .
The slower component (IOSl), however was negative in both transmitted and reflected experiments, thus it can be attributed to the increase of the absorption of the tissue. While different components of the in vitro, blood-free IOS signal were distinguished previously 1,5,20 all of them were attributed to the change of the scattering of the tissue, thus none of them corresponds to the observed increase of the absorption. Cell death could be one of the possible reasons of this phenomena, however, as the electrical responsiveness of the slices haven't changed www.nature.com/scientificreports www.nature.com/scientificreports/ significantly during the recording, we can suppose, that cellular destruction is not necessarily substantial. In order to verify the possible role of the cell death in the IOSl component, its irreversibility should be checked.
A candidate mechanism, showing that absorption plays role in IOS generation, was presented by Mané and Müller 21 : Multispectral analysis of the IOS showed a significant dip in the IOS spectrum at 440 nm wavelength during spreading depression. This absorption line, referred to as Soret band, is a mark of the absorption of porphyrins. The changes of the light absorption by cytochromes appear at this wavelength, thus, reduction of cytorchromes as well as its unpacking from the mitochondria could increase the absorption of the tissue.
The directed causal relationships were examined between the LFP and the IOSh by Sugihara's causality analysis method, the convergent cross mapping. We have found a strong, unidirectional, delayed causal effect from LFP to IOSh with 0.5-1 s delay, without signs of feedback from the IOSh to the LFP. However, it is also showed, that a "shadow peak" in the IOSh → LFP direction appeared at anticausal time delays.
In general, we found, that the delay of the causality peaks significantly affects the interpretation of the results. Assuming no observational delay, two different cases should be distinguished: the peak of the cross-map functions, located on the negative half of the time axis, should be considered as a sign of the real causal effect, since the cause preceded the consequence. However, the peaks on the positive half-axis should be interpreted, as a sign of the delayed prediction, ie. the caused time series follows the effect of the cause faithfully with some delay, thus not only the cause can be reconstructed from the caused time series (this corresponds to the peaks on the negative axis), but the caused time series can be predicted from the cause as well. While the peaks on the negative half of the delay axis are signs of the causality in terms of the Sugihara et al. 12 , the peaks on the positive half-axis correspond more to the predictive causality according to the Wiener-Granger principle 22-24 . In our case, due to the relatively smooth and continuous nature of the IOSh signal, the IOSh → LFP causality peak on the anti causal half-axis become wider, and the "shoulder" of that peak produced relatively high causality values at the negative half of the time lag axis as well. If only the instantaneous CCM were calculated, the IOSh → LFP causality coefficients would have similar values to the LFP → IOSh drive, thus one would erroneously conclude, that a significant feedback exists in the IOSh → LFP direction as well. Observation of the full, time dependent, cross-map function reveals, that the relatively high causality at zero time lag in the IOSh → LFP direction is only a side effect of a peak in the anti-causal regime. We concluded, that the delay dependent calculation of the CCM is important not only to determine the delays of the causal effects, but to find the correct causal directions and strengths as well.
The causality structure was very similar in all the three induction methods of the epileptic activity: electrical stimulation, Mg 2+ -free solution or 4AP induction, although the causal connection strength was stronger in the stimulated and weaker in the presence of 4AP than in Mg 2+ -free solution. The similarity of the causal structures may imply similar underlying mechanisms in all three cases, meaning, that the generation of the IOSh component is independent of the method of induction.
It was demonstrated, that although the cross-correlation functions showed peaks, those peaks did not reflect the actual causal dependency in these cases. Instead, the temporal derivative of the IOSh was correlated with the LFP power at the time delay of the actual causal peak.
During the interdischarge intervals, the IOSh signal decays towards the baseline exponentially, without significant causal influence from the LFP, which implies a linear autonomous dynamics.
To sum up our observations on the dynamics of IOSh, a simple model has been set up to describe the dependency of the IOSh on the LFP power. The model allowed the reconstruction of the IOSh based on the LFP signal. Model fitting showed, that during the slow development of the epileptiform activity, the effective causal connection strength slowly decreased. This decrease could be well approximated by an exponential decay. Quantitatively good reconstruction of IOSh, based on the LFP signal by our model during the 1 h long recordings, supports the results of the causality analysis as well.
The model provided us the possibility to cross-validate the causality analysis, with known causal dependency and effect delays. The causality analysis on the simulated data series resulted in very similar cross-map functions on the model as it was found on the measurements: the peaks on the negative half of the time delay axis precisely marked the direction and the delay of the simulated causal effects, while "shadow peaks" were generated on the positive, anti-causal half of the delay axis. The only slight difference was that in case of simulated data series, the shadow peaks position was always symmetric to the main peak on the delay axis, while in case of the measurements, the "shadow peak" showed larger lags.
Our results on LFP → IOSh unidirectional causality do not exclude the possible presence of causal connections on a slower time scale either in the reverse direction or in the same direction but between the observable epileptic bursts, where we now observe only the autonomous dynamics of the IOSh. Inference of possible interactions on slower time scale now excluded by filtering out the slower component (IOSl) which was necessary to make the IOSh available for the causal analysis.The causal relation between the LFP and the slower, IOSl component could not be determined based on these measurements, because the IOSl does not satisfy the necessary conditions. The application of the CCM method requires, that the system went through multiple times on its attractor. This condition is satisfied for the IOSh component, but not for the IOSl, which exhibits only one sweep through the state space, during these experiments.
The application of CCM for LFP and IOSh demonstrates, that the significantly different inherent speed and sampling rate of the signals, and the necessary downsampling does not preclude the determination of the causal relationships, thus raises the possibility of application of the new method to other signal modalities with different speed, such as fMRI and EEG as well.

Slice preparation.
In vitro experiments were performed on rat cortical slices (13 slices from 13 rats) using standard procedures. Animals were decapitated in deep chloral-hydrate (Hungaropharma, Budapest, Hungary) anesthesia and the brain was quickly removed from the skull. A block was dissected from the somatosensory cortex and 400 μm thick coronal slices were prepared with a vibratome (EMS-4000, Electron Microscopy Sciences, Fort Washington, PA, USA) in ice-cold artificial cerebrospinal fluid (ACSF). Slices were incubated at room temperature for an hour in oxygenated ACSF (pH 7.01-7.12), the composition of which was (in mM): 126 NaCl; 1. Induction and recording of seizure activity -interface conditions. A single slice was placed into an interface-type recording chamber, which was perfused continuously (3 ml/min) with ACSF. In the recording chamber, slices were maintained at 32 ± 0.5 °C in humidified, carbogenated (95%O 2 , 5%CO 2 ) gas atmosphere. A total number of 12 recording sites in 12 slices were included in the final analysis. To test the change in the LFP, glass extracellular electrode (5-10 MΩ) filled with 1 M NaCl was positioned into layer 3 of the somatosensory cortex. For stimulation a bipolar tungsten electrode was placed at the border of the white and gray matter below the recording electrode (Fig. 1). Electrophysiological signals were amplified by an Axoclamp 2B amplifier (Axon Instruments Inc., Union City, CA), filtered and further amplified (0.16-1000 Hz, 1000x) by a Supertech Signal Conditioner (Supertech Ltd., Pécs, Hungary) and digitized by an NI-6023E A/D card (National Instruments, Austin, Texas) for off-line analysis. Intrinsic optical signals (IOS) were recorded also from the beginning of treatment with convulsant solution, simultaneously with local field potential (LFP). The slice was illuminated with unfiltered white light, using a voltage-stabilized cold light source (Fiber-Lite MI-150, Dolan-Jenner, Boxborough, UK). Alteration of reflected light was recorded by a monochrome 12-bit CCD camera (FOculus FO-432B, NET Gmbh, Lerchenberg, Germany) attached to an upright 3-way Olympus SZX-9 (Olympus, Tokyo, Japan) stereomicroscope. Digital images were taken at a resolution of 1024 × 768 pixels at 2 fps sampling rate, and stored in uncompressed format on an x86 based personal computer for off-line analysis Epileptiform discharges were provoked either by high frequency electric stimulation or by convulsant application into the perfusion solution. Evoked seizures developed as afterdischarges, which immediately followed the high frequency electric stimulation (50 Hz, 5 s). Spontaneous epileptiform activity developed after exchanging the perfusion solution in the recording chamber from ACSF to Mg 2+ -free-or 4-aminopiridine (4-AP, 50 μM) containing solution. These types of activity usually appeared in 10-25 min in following the perfusion of the convulsant.
Data acquisition of IOS and LFP was made simultaneously by a custom Matlab (The MathWorks Inc., Natick, MA, USA) based software, which made snapshots of continuous video signal from the CCD camera at a pre-defined sampling rate, based on the timing process of NI-6023E A/D card. For the analysis of the optical changes a region of interests (ROI) was manually defined in the slice beside the site of electrophysiological recording, optical alterations in a small box was determined (Fig. 1).
Recording of seizure activity -submerged conditions. Under submerged conditions the recording of LFP was carried on a MEA USB-1060 INV (Multichannel Systems, Reutlingen, Germany) system. Slices were transferred to a 60 channel 200/30 3D multielectrode array (MEA) chip (Qwane Biosciences, Lausanne, Switzerland), where 30 µm diameter electrodes are positioned in a 8 × 8 grid of 200 µm distances. The slice is positioned directly on the bottom of MEA chip using a net, electrodes with the height of 50-70 Âµm penetrate into the lower surface of slice. Continuous ACSF perfusion was carried out at a 5-6 ml/min rate, temperature is maintained at 33 ± 1 °C using TC02 Temperature Controller (Multichannel Systems, Reutlingen, Germany) with PH01 perfusion heater (Multichannel Systems, Reutlingen, Germany) and the built-in heating pad of MEA 1060 Amplifier (Multichannel Systems, Reutlingen, Germany). The slice was trans illuminated with unfiltered white light, using a Rebel Star LED (Ledium Ltd, Szeged, Hungary) light source, placed under the MEA chip. Alteration of transmitted light was recorded by a 10-bit CCD camera (Qimaging Micropublisher 3.3, Qimaging, Surrey, Canada) attached to an upright 3-way Nikon SMZ 800 (Nikon Instruments Europe BV, Amsterdam, Netherlands) stereomicroscope. Digital images were taken at a resolution of 1024 × 768 pixels at 1 fps sampling rate, and stored in uncompressed format on an x86 based personal computer for off-line analysis Data acquisition of IOS and LFP was also made simultaneously in the case of submerged slice experiments. The STG4002 stimulator (Multichannel Systems, Reutlingen, Germany) provided continuous triggering for a custom Matlab (The MathWorks Inc., Natick, MA, USA) based framegrabber software during data acquisition made by MEA 1060 amplifier. Epileptic activity was provoked by Mg 2+ -free perfusion solution application.
The intrinsic optical signals. Image series of interface and submerged slice experiments were analyzed with the same method, using a custom Mathworks Matlab based software. Since the first 2-3 min of video data acquisition was made in normal ACSF, the first 10 frames of each image series were averaged and served as control F cont , and then it was subtracted from each subsequent images according to the following formula: where F cont is the luminance of control image (the average of first 10 frames) at a given pixel, F t is the luminance of each subsequent experimental image, ΔF is the change of reflectance or transmittance, depending on experimental conditions. F cont in the denominator serves to normalize the data across regions differing in luminance. The images were converted to 256-scale pseudocolour images for further off-line analysis. For causality analysis a region of interests (ROI) was defined in the somatosensory cortical region of the slice around the site of electrophysiological recording. The mean IOS value was calculated for the ROI. This data was used to characterize the IOS changes of the slice during 1 h seizure induction and correlate with single-channel field potential recordings.
www.nature.com/scientificreports www.nature.com/scientificreports/ Ye et al. 13 extended the original CCM method by introducing time delay between the cause and consequence and examining the causal relationships at many time lags.
Theoretically, if the dimension of the dynamics is D, than at most an E = 2D + 1 dimensional embedding space is needed to reach the full topological equivalence between the embedded and the original attractors 25 . However, E = D + 1 embedding dimension is enough, to turn the false self-cross sections of the manifold into zero-volume 26 . In practice, we found that higher embedding dimension and higher τ embedding delay resulted in increasing causality-peak for a while, but widening plateau on the top of the peak after a certain point. Thus, we tried to find parameter set, which is optimal for the inference of the effect's temporal delay, which results in the highest peak without plateau on it. As a result, the embedding dimension was chosen to 3 and the time delay for the embedding to 1 time step (0.5 s) for the analysis.

Data Availability
The datasets generated during and analyzed during the current study are available from the corresponding author on request.