Dopamine receptor antagonists effects on low-dimensional attractors of local field potentials in optogenetic mice

Summary

The goal of this study was to investigate the effects of acute cocaine injection or dopamine (DA) receptor antagonists on the medial prefrontal cortex (mPFC) gamma oscillations and their relationship to short term neuroadaptation that may mediate addiction. For this purpose, optogenetically evoked local field potentials (LFPs) in response to a brief 10 ms laser light pulse were recorded from 17 mice. D1-like receptor antagonist SCH 23390 or D2-like receptor antagonist sulpiride, or both, were administered either before or after cocaine. A Euclidian distance-based dendrogram classifier separated the 100 trials for each animal in disjoint clusters. When baseline and DA receptor antagonists trials were combined in a single trial, a minimum of 20% overlap occurred in some dendrogram clusters, which suggests a possible common, invariant, dynamic mechanism shared by both baseline and DA receptor antagonists data. The delay-embedding method of neural activity reconstruction was performed using the correlation time and mutual information to determine the lag/correlation time of LFPs and false nearest neighbors to determine the embedding dimension. We found that DA receptor antagonists applied before cocaine cancels out the effect of cocaine and leaves the lag time distributions at baseline values. On the other hand, cocaine applied after DA receptor antagonists shifts the lag time distributions to longer durations, i.e. increase the correlation time of LFPs. Fourier analysis showed that a reasonable accurate decomposition of the LFP data can be obtained with a relatively small (less than ten) Fourier coefficients.

Keywords:

Prefrontal cortex – Neural networks – Lasers – Cocaine – Pyramidal cells – Fourier analysis – Interneurons – optogenetics

Introduction

Use of psychostimulants, such as cocaine, is a serious health problem and opens the door to neurobiological changes in limbic and cortical circuits that engage cognitive and emotive processing. We have just began to understand the cellular adaptations that occur in the cortex following a single exposure to cocaine and their contribution to the continuous and further use of drugs of abuse. The behavioral consequences of first time cocaine use are varied and appear to be somewhat contradictory. First time cocaine users often report feeling a sharpening of the senses [1], and anecdotal information suggest that acute cocaine increases attention. Indeed, individuals with ADHD will sometimes self-medicate with cocaine [2]. Contrastingly, Jentsch and colleagues [3] have shown that acute cocaine administration impairs performance on a reversal learning task, and several studies have reported compromised performance during repeated acquisition tasks in monkeys [4, 5]. Additionally, imaging studies in humans have shown that acute cocaine administration induces prominent activation of the prefrontal cortex, primarily in the dorsolateral regions [6]. Furthermore, acute cocaine administration has been linked to poor impulse control [3, 7, 8]. Therefore, it seems that first time cocaine use may give users a sense of enhanced awareness, while cognitive performance is diminished. Cocaine abuse is a public health challenge for which we are still seeking more effective treatments [9, 10]. Additionally, patients relapse in greater numbers increasing the social and economic cost of this disease and the burden on medical and welfare systems. Thus, understanding the cellular and network process underlying the effects of cocaine in the prefrontal cortex is a critical necessity.

Neuronal oscillations are thought to be general and fundamental mechanisms for enabling coordinated activity during normal brain functions [11–13] and it is within this framework that the function of gamma oscillations have recently acquired importance. Gamma oscillations in the cortex involve the reciprocal interaction between interneurons, mainly PV+ fast spiking interneurons (FS PV+) and principal cells [14]. Gamma oscillations appear to be a critical mechanism underlying the cognitive and behavioral function of mPFC. It is therefore highly likely that gamma oscillations in mPFC would be altered by cocaine administration. In this study we investigated changes in gamma oscillations following an acute administration of cocaine and the effects of selective dopamine (DA) receptor antagonist on the gamma band oscillations.

Parvalbumin (PV) interneurons and pyramidal cells in the mPFC

The vast majority (over 70% according to [15, 16]) of excitatory neural population is made of pyramidal cells. Pyramidal cells span cortical layers 2-6 [16] and are thought to be involved in both maintaining the working memory [17] and the UP states of persistent network activity [18, 19]. The main modulators of the pyramidal cells are interneurons, of which PV+ fast-spiking interneurons represent the majority [20–22]. PV+ neurons coordinate the output of the local minicolumns [23, 24], maintain and/or modulate both the 15–30 Hz beta [25–27] and the 25–40 Hz gamma [28–30] rhythms of the brain, and facilitate information processing [31–33]. PV+ interneurons are known to both modulate the firing rate of pyramidal cells through their dendritic projections [34, 35] and to block action potentials through their axonal projections [34–38]. Abnormal activity of PV+ inerneurons have been linked to autism [39–42], schizophrenia [32, 43–45], sensory hypersensitvity and neural hyper-excitability [46–49].

Optogenetic studies of mPFC

Optogenetics refers to the use of optical stimuli to elicit neural responses from neurons whose ionic channels were genetically modified to respond to light stimuli. Optical fibers transmit the light stimulus from the external laser to the brain region of the mouse. For optogenetics, the optical fiber is integrated with electrophysiological probes and form a device called an optrode [50]. In basic science studies, optogentics has been used for investigating neural plasticity mechanisms [51–53], understanding information processing by neural networks [29, 30, 54], and testing hippocampal memory formation hypotheses [55, 56]. In behavioral studies, optogentics proved invaluable in studying feeding [57–60], fear conditioning [61, 62], and aggression [63]. The ultimate goal of these studies is to find optogenetics-based solutions to neural disorders, such as anxiety and depression [64–66], and ameliorate neurodegenerative conditions, such as Parkinson disease [67, 68], and epilepsy [69–72]. Finally, another promising development includes the use of optogentics for the restoration of visual functions in blind animals [73–75] and the design of neural interface [76].

Dopamine (DA) receptor

The widely-accepted paradigm regarding the mechanisms of cocaine and amphetamine effect is through increasing synaptic levels of DA, mainly through DA transporter’s inhibition [77–80]. As a result, DA receptors are exposed to an elevated and sustained levels of endogenous DA, which amounts to an indirect DA agonist effect [81]. All drugs of abuse produce increased levels of DA [82–85]. Additionally, the relapse to cocaine-seeking behavior implicates the mesolimbic DA system [86–89] through DA receptor’s potentiation [86, 90].

There are several kinds of DA receptors, which are studied with the purpose of finding effective addiction treatments [91–93]. Five DA receptor subtypes, called D1-D5, have been identified [94–97], although they appear to form two classes of pharmacologically and biochemically similar receptors, i.e. D1-like receptors, which refer to D1/D5, and D2-like receptors, which refer to D2/D3/D4, respectively. Both D1- and D2-like receptors have both been implicated in the discriminative and reinforcing effects of cocaine [78, 98, 99], although they seem to be somewhat selective in their response to cocaine [100, 101].

We previously carried out a similar computational study on a subset of the same animals under control conditions [102] and following an acute cocaine injection [103]. The present study expands our previous investigation of cocaine effects on mPFC [102, 103] in which we systemically injected DA receptor antagonists & cocaine. As in our previous studies, we corrected for the light-induced phase resetting using both the autocorrelation [102, 104] and Hilbert’s transform [105]. Delay-embedding method has been used for estimating the delay (lag) time and embedding dimension of the LFPs [102, 106, 107]. We also used spectral analysis to model the LFPs and found that a significant fraction of the data trend can be captured with as few as ten Fourier coefficients.

Materials and methods

A detailed description of the procedures can be found in the previous two papers of this series [102, 103]. Mice Cre-PV+ adult males were mounted in an steretoaxic apparatus, and a small hole drilled over the PFC to insert the optrodes. Temperature was maintained at 35 C via a thermal blanket and isoflourane anesthesia was maintained using an isoflourane anesthesia system. At the end of the recordings animals were deeply anesthetized and the brain extracted to confirm anatomical localization of optrodes. All procedures were done in accordance to the National Institute of Health guidelines as approved by the Medical University of South Carolina Institutional Animal Care and Use Committee.

Experimental protocol

The response of medial prefrontal cortex (mPFC) to a brief 10 ms optical stimulation was measured from 2-s long LFP recordings that were repeated 100 times for each of the 17 animals in this study. The fast-spiking interneurons were activated by laser stimulation. This study used the same mice model together with optogenetics and in vivo electrophysiology described in detail in [108]. The viral vector (AAV2/5. EF1a. DIO. hChR2(H134R)—EYFP. WPRE. hGH, Penn Vector Core, University of Pennsylvania) has been injected to the mPFC of male mice at least 3 weeks before electrophysiological recordings (B6; 129P2 − Pval^{btm1(Cre)Arbr/J} Jackson Laboratory (Bar Harbor, ME, USA) [108]. The extracellular LFPs were amplified (Grass Technologies, West Warwick, RI, USA), digitized at 10 kHz (1401plus data acquisition system), and analyzed offline. The signal was filtered (Quest Scientific Inc., Canada) and band-pass online between 0.1 and 130 Hz to obtain the LFPs. The light stimulation was provided by a 473 nm laser (DPSS Laser System, OEM Laser Systems Inc, East Lansing, MI, USA).

The animals received acute systemic injection of the D1 receptor antagonist SCH 23390, 1.0 mg/kg, i.p. or the D2 receptor antagonist sulpiride, 15 mg/kg, i.p. and optical stimulation.

The novelty of this study compared to others involving DA receptor antagonists is the use of delay embedding for data visualization and spectral analysis for data modeling. The novelty of this study compared to the previous two [102, 103] is that we targeted DA receptors with specific receptor antagonists. When the DA receptor antagonist was injected systemically after the baseline and prior to cocaine the name of the trial is xxxbyy where xxx is either sch for SCH 23390 or sulp for sulpiride, the letter b stands for before cocaine and yy is the mouse number. For example, schb14 means LFP recordings from mouse # 14 with SCH 23390 D1 receptor antagonist applied before cocaine. Similarly, the name xxxayy means the receptor antagonist xxx was applied after cocaine. Finally, in other trials we applied both DA receptor antagonists together and the name of the trial is both followed by either a (after cocaine) or b (before cocaine) and the mouse number (see Table 1). For example, bothb24 means LFP recordings from mouse # 24 with SCH 23390 & sulpiride applied before cocaine.

In the following, we briefly summarize the main methods used here for data analysis, visualization, and modeling.

Cluster analysis with dendrograms

A dendrogram is a diagram that shows the hierarchical relationship between objects. The dendrograms are used for allocating objects to clusters. For example, a scattered plot in a two-dimensional Euclidian space may look like Fig 1a. The (x, y) pair could be single measurements data or, as in our study, they could represent an entire trial. In our example, the dendrogram computed the “distance” between the scattered points shown in Fig 1a and then ordered them according to the distance, which is shown along the vertical axis of the dendrogram in Fig 1b. In the example above, we can immediately notice from Fig 1a that data points #4 and #5 are the closest across all ten, and therefore they will be shown in the dendrogram of Fig 1b as the most similar, i.e. the height of the link that joins them together is the smallest. A cursory inspection of Fig 1a suggests that the next closest data points are #1 and #2, which translates into the next most similar objects in the dendrogram shown in Fig 1b. In a dendrogram, the height of the dendrogram indicates the order in which the clusters were joined. Each joining (fusion) of two clusters is represented on the dendrogram by the splitting of a vertical line into two vertical lines. The vertical position of the split, shown by the short horizontal bar called clade, gives the distance (dissimilarity) between the two clusters. The terminal end of each clade is called a leaf.

Observations are allocated to clusters by drawing a horizontal line through the dendrogram. In our example, below the dashed line marked (a) in Fig 1b there are two clusters, which are separated by the corresponding dashed line in Fig 1a. If the required maximum distance is lowered (see dashed-dotted line marked (b) in Fig 1b), then we identify three clusters, which correspond to separated data clusters in Fig 1a. By further decreasing the required maximum distance within a cluster (see the double dashed line marked (c) in Fig 1b), then five clusters emerge.

Although in our example, we simply used the Euclidian distance between the points shown in Fig 1a, other distance measures could be used for hierarchical organization of data. Among many possible distance measures, the correlation coefficient [109] and the Euclidian distance [110] are most common. For example, a correlation-based distance is distance = 1 − c, where c is the correlation coefficient (c).

Delay-embedding method

In this paper series (see also [102, 103]), we used the delay-embedding method of nonlinear dynamics to estimate the number of degrees of freedom of the steady activity of mPFC neural network. The challenge is to extract, or reconstruct, from the one-dimensional data (time series) the actual multidimensional attractors of neural activity.

Theoretically, the first derivative of continuous functions, or equivalently the first order difference of discrete data sets, is orthogonal to the original time series. Therefore, one could use the original time series and its first order difference as the two “natural” dimensions (axes) to unfold the dynamics of the system. If the phase space trajectories cross each other, then we could use the second derivative as the third orthogonal direction in the phase space and try to unfold the attractor again. The procedure could be repeated until all self-crossings of the phases space trajectory are eliminated, which is the final unfolded attractor of neural dynamics [111]. Takens’ theorem proved that for regularly sampled and noise free time series, the attractor reconstructed using the above-described differential embedding is diffeomorphic to the original (unknown) attractor [112, 113].

A more robust to noise approach is to use the delay-embedding method [112, 114], which takes a time series x_i = x(iΔt) with i = 1, 2, …,N where N is the number of data points and Δt is the uniform sampling time, and expands it into a d−dimensional vector by temporally shifting the data [112, 114–116]:

The lag time

The challenge of selecting the “right” delay time is to avoid both the redundancy, i.e. high correlation among the data points because of a too small delay time, and completely de-correlating the data points by a too large delay time [117]. As in the previous studies [102, 103], we used both the autocorrelation [117–122] and the average mutual information (AMI) [123–125] for estimating the lag time τ. The first zero crossing of the autocorrelation ensures that the lag time completely de-correlated the original and the time-shifted time series [126]. While convenient, the autocorrelation method only eliminates the linear correlation. As a result, we also estimating the lag time using the first minimum of the nonlinear autocorrelation function called Average Mutual Information (AMI) [123, 125].

In order to allow a direct comparison of our current results on the effect of DA receptor antagonist combined with cocaine on LFP recordings against similar results on the same mice for baseline recordings (see [102]) and cocaine (see [103]), we only analyzed the lag times obtained from the autocorrelation function method. However, for completeness, we also briefly summarize the statistics of the lag times obtained from AMI.

The embedding dimension

As before [102, 103], we used the false nearest neighbors (FNN) algorithm [124, 127, 128] to estimate the embedding dimension. Intuitively, a high-dimensional phase space trajectories projected onto a too low dimensional embedding space has self-crossing points, or false nearest neighbors (FNN). Such FNN could be eliminated by progressively increasing the dimension of the embedding space until reaching a reasonably low percentage of FNN [127].

Without repeating all the details in [102, 103], we used the lag times computed in the previous section for each condition to estimate the embedding dimension with variable distance ratios, f, between 2 and 20. The distance ratio f estimates the factor by which the distance between two points increased when increasing the embedding dimension by one unit [126, 127, 129]. Usually, small f overestimates the percentage of FNN, whereas a too large f values give a large number of false positives. In our study, for f > 10, the percentage of FNN drops below 1% for an embedding dimension d_E = 3 (not shown), whereas for f > 12 the percentage of FNN dropped below 10⁻⁶% at d_E = 3.

As in the previous studies [102, 103], we further removed spurious temporal correlations [130, 131] with Theiler windows from 100 to 8000 sampling times [131, 132]. For all tested Theiler windows, the percentage of FNN dropped below 0.1% for an embedding dimension of d_E = 3 and f > 10.

Fréchet distance method

Fréchet distance [133] relates to the “man walking a dog on a leash” paradigm, i.e. find the shortest leash needed such that the man walks along one curve and the dog walks along the other [134]. Both the man and the dog can adjust their speeds but are not allowed to move backwards [135]. Fréchet distance is a topological measure of the similarity between two curves [134], i.e. it does not consider any temporal constrains regarding the location of figurative points on the two curves. We computed Fréchet distance using the algorithm in [136].

Fourier method

We used Fourier analysis to determine the spectral composition of each trial. Any periodic signal x(t) of period T could be represented as an infinite sum of harmonics [137]:

Results

As in the previous papers of this series [102, 103], the response of the local mPFC neural network to a brief 10 ms light pulse was recorded for 2 s and the procedure repeated 100 times for each animal. We discarded the transient response of the neural network by removing the first 0.5 s of each trial. The transient response contains information regarding the phase resetting of the entire local neural network, which has been successfully applied to epilepsy [139, 140] or Parkinson’s [141, 142] studies. However, here we focused on the last 1.5 s of steady activity of each trial with the purpose of identifying stable and repeatable patterns of neural activity that may connect the control [102], the pure cocaine [103], and this DA receptor antagonists study.

Phase resetting of LFPs

A brief 10 ms light pulse alters both the phase and the amplitude of the ongoing oscillatory activity of the local mPFC [104, 106, 143, 144]. Such transient resettings dissipate after a few cycles and eventually the neural network returns to its unperturbed steady activity [105]. We corrected the trials for the transient phase resetting induced by light stimuli by circularly shifting each LFP trace with respect to an arbitrary “reference” trial until maximizing the correlation coefficient [105, 144]. Phase resetting correction led to a significant increase in the coefficient of correlation among trials as seen in Fig 2a1 for the case of scha16 and in Fig 2a2 for all mice and conditions. The correlation coefficients before circularly shifting the data to account for the transient LFP phase resetting has very low values (see the solid black trace in Fig 2a1 and solid black squares in Fig 2a2). For each condition, e.g. scha, schb, supla, sulpb, botha, bothb, and each animal #, the averages of all correlation coefficients (such as the one shown in panel a1) indicate a very low initial correlation (solid black squares in Fig 2a2) due to the fact that the incoming light stimulus finds the neural network at random phases. By circularly shifting the trials to maximize the correlation coefficient we corrected for the random phase of light stimulus relative to the ongoing rhythmic activity of the network (solid red circles in Fig 2a2). Although for the same animal and the same condition the correlation coefficient differs from one pair of trials to the other as shown in Fig 2a1, we averaged the correlation coefficients over all 100 trials and only showed in Fig 2a2 the average and standard deviation bar. We further averaged all data for the same condition to gain a better understanding of the average increase of correlation due to phase shifting.

For example, the three scha16, scha24, and scha25 data points in Fig 2a2 (already averaged over all trials) were considered as a single scha condition in Table 2. While the improvement in the correlation coefficient is qualitatively apparent from Fig 2, the compact numerical summary in Table 2 shows that in most cases the circular shifting improvement exceeded one order of magnitude.

We also computed the sum of squared differences (SSD) between the “reference” trial and all the other trials for each animal and each condition (see one example in Fig 2b1 for scha16). As expected, circularly shifting the trials to maximize the correlation coefficient also reduces the SSD among trials (Fig 2b1), overall per animal and condition (Fig 2b2), and also when aggregated in a single average per condition as in Table 2.

Dendrograms of phase shifted LFPs

A dendrogram is a visual representation of “relationships” among trials. The horizontal axis of the dendrogram in Fig 3a1–3b1 shows the trial index. However, the trials are not listed in their recording order but rather based on their membership to the respective cluster. The vertical axis of the dendrogram is the distance between the “leaves”, i.e. trials, of the dendrogram tree (see Fig 3a1–3b1).

We used dendrograms for two purposes: (1) to qualitatively check that the circular shifting of trials for the purpose of maximizing the correlation coefficient also reduced the distance among trials in a cluster, and (2) to quantitatively measure the amount of overlap between the baseline and the DA receptor antagonists trials.

Since we established that the correlation coefficients increase by circularly shifting the trials (see Fig 2a1 and 2a2) then it follows that the corresponding dendrogram distances would decrease (see Fig 3a1–3b1). In particular, it is expected that the distances in the circularly shifted dendrogram (Fig 3b1) should be smaller compared to the original data (Fig 3a1). The reason we used the dendrogram of the circularly shifted data is because it better separated trials in disjoint classes as we can see from two examples of LFP traces shown Fig 3a2 and 3b2. Indeed, before circularly shifting the data, even trials as dissimilar-looking as those shown in panels a2 and b2 could be classified as belonging to the same cluster in Fig 3a1. This is because a certain amount of correlation, no matter how small, would still exist between two trials recorded from the same system even if they have a random phase difference between them due to light stimulus resetting of neural activity. By circularly shifting the trials with respect to one (reference) trial to maximize the correlation coefficient, we eliminated random overlaps among trials which resulted in small correlations shown in Fig 2a1 and 2a2 and only retained for further analysis correlations that are significant, i.e. one order of magnitude larger than before circularly shifting the data.

The Euclidian distance that measures the similarities among the “leaf” nodes of the dendrogram (see Fig 3a1 and 3b1) has been significantly reduced by phase reseting correction. For example, for the scha16 experiment the distance threshold for breaking the 100 trials of into six clusters was reduced from 62 (arbitrary) units (see the horizontal dashed line in Fig 3a1) to 42 units (see the horizontal dashed line in Fig 3b1).

In Fig 3a2–3b2 we plotted two traces from the first two largest clusters of the 12 shown in Fig 3b1. There are clear differences between the traces, e.g. the largest cluster in Fig 3b1 collects all LFP with small amplitude and high frequency oscillations (Fig 3a2), whereas the second largest cluster collects large amplitude LFPs that show a graded, almost linear, amplitude decay over time (Fig 3b2).

Dendrogram-based quantitative test

We also used the dendrogram-based clustering for measuring the percentage overlap between baseline and DA receptor antagonists trials. For each animal, we concatenated the 100 baseline trials (trial index from 1 to 100) with the 100 trials after DA receptor antagonist treatment (trial index 101 to 200) in a single 200-trial file. If the neural activity patterns for baseline and DA receptor antagonists were totally distinct, then we would expect that they separate in disjunct clusters. As a result, a perfect dendrogram classifier would correctly separates the 200 mixed baseline and dopamine antagonists trials in disjoints (non-overlapping) clusters such that each of them only contains one type of trials, either only baseline or only DA receptor antagonist trials. On the other hand, if half of the trials in a given cluster are from baseline and the other half are from DA receptor antagonists then the overlap is maximum possible and the Euclidian distance-based dendrogram cannot discriminate between those trials.

For each cluster, we computed the percentage of mixing between baseline and DA receptor antagonists trials by using the formula [102, 103]:

On the low side, the 20% average mix between trials could indicate the lower bound of the correlation-based dendrogram sensitivity, which classifies some trials as “similar”, regardless how small is the correlation among them. The fact that even after circularly shifting the data we still find 20% misclassification between baseline and DA antagonist trials is an indication that the correlation-based clustering method is not perfect. The 50% average mix suggests that for the respective mice (see Fig 4c) the correlation-based dendrogram classifier is not effective in distinguishing between the baseline and DA receptor antagonist. For this purpose, we further investigated the statistics of lag times using the Fréchet distance among trials to sort out true similarities (see below).

Lag time statistics

The autocorrelation-based lag time distributions for all six conditions listed in Table 1 are shown in Fig 5. All weighted averages of lag times were presented with four significant figures because the sampling rate was 10 kHz, i.e. a sampling time of 0.0001 s.

SCH 23390 after cocaine (Fig 5a1) shows a clear bimodal distribution of lag times with one distinct peak at very low lag times (around 0.1s) and another peak at much longer lag times (about 0.45 s). For every pharmacological manipulation, such as sch, sulp, or both from Table 1, there is always a baseline recording for the respective mouse [102]. Such baseline recordings were called bcxx where xx is the mouse number and the corresponding distributions of lag times are shown in (Fig 6). In order to quantify the effect of DA receptor antagonists on cocaine administration, we could directly compare the distributions of lag times against the baseline distributions. For example, Fig 5a1 (SCH 23390 condition mice # 16 and 24) directly compares against Fig 6a1 (baseline recordings mice # 16 and 24). While the distribution of lag times under scha treatment, i.e. cocaine after SCH 23390, preserves the original bimodal distribution of the baseline it has a noticeable drift towards longer lag times. This means that the correlation time of the scha data increased compared to baseline.

Another direct comparison is with the cocaine data [103] shown in Fig 7. As we previously reported, cocaine generally shifts the lag times towards shorter durations compared to baseline values [103]. If we compare the peak around the lag time of 0.1 s in the baseline Fig 6a1 against the cocaine Fig 7a1, then not only that the peak narrowed under cocaine but it is also much stronger compared to any other lag times. More quantitatively, the weighted average of all lag times for the baseline was 0.2072 s, 0.1937 s after cocaine, and 0.3012 s after scha (see Table 3). In this study, we expanded our previous analysis [102, 103] from six to 17 animals and noticed that sometimes the cocaine increases the weighted average of lag times compared to baseline. However, to put this apparently contradictory result in statistical perspective, we notice that the largest increase of 35% is only observed in one animal, i.e. bc23. The other increase of 17% compared to the baseline is also due to only one animal bc10 that skews the entire average of the group. In all the other 15 animal the cocaine shifts the weighted average lag times to shorter durations compared to the baseline as previously reported in [103].

This confirms the observation that generally cocaine alone shifts the lag times towards smaller durations compared to baseline [103], i.e. decreases the correlation time of the data, whereas SCH 23390 applied after cocaine significantly shifts the lag times towards longer durations. In this study, we used the weighted average of the lag times to put a number on a complex and complicated distribution of lag times. In our previous studies, we fitted the lag time distributions with log-normal functions and extracting the peak of the fitting function (see [103]). The main reason for such a change in our approach is that a log-normal fit, while appropriate for the six mice considered in [103], it is not appropriate for all distributions in this 17 mice study. Indeed, some distributions are clearly bimodal, such as the one shown in Fig 5a1, which would be a very poor log-normal fit.

SCH 23390 before cocaine lag time distributions (Fig 5b1) are unimodal with a weighted average of lag times around 0.26 s. The range of lag times is similar for both SCH 23390 conditions and is shorter than 0.6 s. The unimodal shape reflects the original baseline unimodal shape of the lag time distribution in Fig 6b1. The weighted average of lag times for baseline was 0.2594 s (see the first column in Table 3), which shifted to 0.2635 s under cocaine (see Fig 7b1), which remained unchanged when DA receptor antagonists were applied before cocaine (see the second column in Table 3). The shift in the lag time when cocaine follows the DA receptor antagonists is shown in the third column of Table 3. When SCH 23390 was applied before cocaine (second column in Table 3) it facilitated the shift of the lag times back towards the baseline values (first column in Table 3) and canceled out the effect observed in cocaine trials (third column in Table 3). When SCH 23390 was applied after cocaine (fourth column in Table 3) it facilitated the occurrence of longer lag times compared to baseline.

Sulpiride after cocaine (Fig 5a2) also shows a bimodal distribution with one peak around 0.25 s and the other, significant but less prominent, around 0.5 s. The bimodal distribution is more prominent in the baseline histograms shown in Fig 6a2 and also somewhat preserved under cocaine (see Fig 7a2). The weighted average of lag times for the baseline was 0.2275 s, remained unchanged under cocaine (third column in Table 3), and then shifted to a longer duration of 0.2864 s when sulpiride was applied after cocaine (the fourth column in Table 3). We noticed that when either DA receptor antagonists sulpiride or SCH 23390 are applied after cocaine (see the fourth column in Table 3), they increase the lag time compared to baseline (see the first column in Table 3).

Sulpiride before cocaine trials (Fig 5b2) have almost unimodal distributions of lag times with a peak around 0.2 s. The range of both distributions is a bit over 0.6 s. The weighted average of lag times was 0.2439 s for the baseline (Fig 6b2), shifted to 0.2850 s under cocaine (Fig 7b2), and returned to the baseline value of 0.2439 s when sulpiride was applied before cocaine (Fig 5b2).

Both SCH 23390 & Sulpiride after cocaine is unimodal (Fig 5a3), although we only retained one good data set. The weighted average of lag times for the baseline was 0.1790 s (see Fig 6a3), 0.2423 s after cocaine (see Fig 7a3), and 0.1790 s for both sulpiride and SCH 23390 after cocaine (see Fig 5a3).

Both SCH 23390 & Sulpiride before cocaine clearly show unimodal lag time distributions (see Fig 5b3). As opposed to SCH 23390 or sulpiride alone, the lag times in this case concentrate around very short values (around 0.12 s) and the distribution has a long tail of significant contributions up to 0.6 s. The weighted average of lag times was 0.2548 s for baseline (see Fig 6b3), 0.2194 s for cocaine (see Fig 7b3), and 0.2548 s when sulpiride and SCH 23390 were applied together before cocaine (see Fig 5a3).

To summarize, regardless of the effect of cocaine when applied alone (the third column in Table 3), the combination DA receptor antagonists followed by cocaine (the second column of Table 3) always leaves the lag time at the baseline values. On the other hand, when DA receptor antagonists is applied after cocaine (the fourth column in Table 3) it increases the lag time compared to the baseline. There is one exception to this observation (see the last line in Table 3), but it could be an outlier since there was only one set of data available for this particular condition.

Below we briefly summarize the results on the lag times obtained from AMI in Fig 8 and Table 4. When comparing the lag time obtained with the autocorrelation method (see Fig 5) against AMI lag times (see Fig 8) one obvious difference is that AMI lag times are about an order of magnitude shorter.

One benefit of the AMI-based lag times analysis is the identification of a few faulty couplings between the optrode and the tissue. Lag times of the order of 1-5 ms could be realistic when studying small neural networks, since they are of the same order of magnitude as the durations of individual action potentials. However, such short lags are unrealistic when measuring LFPs that average processes over hundreds and possibly thousands of individual action potentials. We noticed a few unusually short lag times (between 2 ms and 5 ms) while most of the lag times are above 0.01 s. As shown in Fig 8, such examples are mouse 17 (about 30% of the lag times below 5 ms), mouse 22 (over 85% of the lag times below 5 ms), and mouse 21 (over 40% of lag times below 5 ms). A closer inspection of all these particular trials showed that the LFPs have very low peak-to-peak amplitude (below 0.04) and they were covered by noise, probably due to faulty coupling with the optrode.

One immediate conclusion from comparing autocorrelation-based lag times from Table 3 against the AMI-based lag times from Table 4 is that the standard deviation of AMI-based lag times is smaller in most cases. This suggests a narrower distribution of lag times for AMI-based method. Another conclusion is that the correlation-based lag time (see Fig 5) is larger than the AMI-based lag time (see Fig 8) by an order of magnitude. With such a wide spread of lag times, the concern is that the attractors may not unfold properly. We found that attractors unfold both with the AMI-based and autocorrelation-based lag times.

Kolmogorov-Smirnov (KS) test

We are both interested in determining the exact values of the delay (or lag) time and also check if the lag time distributions differ from one condition to another. This is the second statistical test, in addition to the previously performed dendrogram clustering of similar trials, aimed at identifying statistically significant variability among trials for the purpose of subsequent mathematical modeling. We are interested in checking if baseline trials which previously could not be distinguished from DA receptor antagonist trials using dendrogram classifiers are in any way distinct from each other with respect to lag time distributions.

Although the plots of the distributions in Figs 5–7 offer a good qualitative comparison among them and the wighted average of lag times (see Table 3) gave a quantitative description of the expected lag, we also performed a direct comparison of lag time distributions. Multiple conventional statistical significance test methods can be used, such as two-tailed t-test, two-sample Kolmogorov-Smirnov test, Wilcoxon rank sum test, analysis of variance, and Kruskal-Wallis test.

Kolmogorov-Smirnov test [145] was used here for deciding if a sample comes from a population with a specific distribution. If the data come from the same distribution at the default 5% significance level, then we represented that result with a solid black square in Fig 9a1–9c1. The probability of observing a test statistic as extreme as, or more extreme than, the observed value under the null hypothesis is shown on a natural logarithmic scale in Fig 9a2–9c2. Statistical significance is often referred to as the probability value, or p-value. A small p-value means that the data are unlikely under some null hypothesis. Usually, the null hypothesis is rejected if p < 0.05. We also computed the maximum difference between empirical distribution functions and plotted them in Fig 9a3–9c3.

To understand how much of the pharmacological manipulation changes the distributions of lag times, we first run the KS test on the baseline data (Fig 9a1). Any off-diagonal black square indicates that the respective baseline distributions are from the same empirical distribution. For example, the lag times for bc8, bc9, bc11, bc15, bc20, and bc22 seem to all follow the same distribution.

The next comparison is among the distributions of the baseline versus DA receptor antagonists trials (see Fig 9b1). For example, the mouse #8 has been treated with sulpiride before cocaine (sulpb) and the distribution of delay times for suplb8 is similar to bc8, bc9, bc11, bc15, bc20, and bc22. This means that for this mouse the sulpiride injection did not change the lag time distribution and, therefore, the black squares remained in the same place as in the baseline case shown in Fig 9a1.

For other mice, such as #18 it turned out that the injection of sulpiride after the cocaine (sulpa) changed the distribution of lag times and made it more similar to bc9, bc11, bc12, bc15, bc19, bc20, bc21, and bc22 (see the solid red squares in Fig 9b1). At the same time, for mouse # 18 the sulpa distribution is definitely different than its baseline (bc18) and also different from bc23 (see the solid blue squares in Fig 9b1). Notice that in the baseline case, bc18 and bc23 lag time distributions are of the same type. This suggests that the combination cocaine followed by sulpiride is affecting the dynamic reconstruction of the attractor through the correlation lag time of data.

When compared against the baseline lag time distributions, the changes seem to occur in scha16, sulpa17, sulpa18, sulpa19, and scha24 recordings (solid red squares in Fig 9b1).

Finally, when computing the KS test across all DA receptor antagonist trials (see Fig 9c1) we notice an increase in the fraction of similar distributions compared to the baseline. For example, mouse #24 had similar distributions with mice # 17, 20, and 22 (see solid black squares in Fig 9a1). However, scha treatment made the lag time distributions for mouse #24 similar to those of mice # 12, 16, 18, 19, and 21 (see the solid red squares in Fig 9c1). Notice that before treatment the distribution of lag times for #24 was different from the baseline for # 12, 16, 18, 19, and 21. We also noticed that scha treatment removed similarities with the baseline # 17, 20, and 22, which is shown with green solid squares in Fig 9c1.

Reconstructed neural activity under DA receptor antagonists

The lag time

The distribution of all autocorrelation-based lag times are shown in Fig 5 and the AMI-based lag time values (not shown) were within 5% of those obtained with the autocorrelation method. As before [102, 103], we used the Tisean function autocor to compute the autocorrelation and mutual to compute the AMI.

The embedding dimension

The lag time was used for computing the embedding dimension of the data over distance ratios, f, between 2 and 20 and Theiler windows from 100 to 8000 sampling times (see [102, 103] for explicit Tisean function calls). We found that the percentage of FNN dropped below 0.1% for an embedding dimension of d_E = 3 and f > 10.

The phase space attractor of each trial was obtained but we only show two representative examples from each dendrogram-based cluster (see Fig 10). Although the attractors for different clusters look different, they are topologically equivalent (see [102, 103]), i.e. by changing the delay time one phase space trajectory could be morphed into the other. The attractors in Fig 10a1 (scha), a2 (sulpa), and a3 (botha) show an “8”-shape, which is topologically equivalent with a closed elliptic attractor as shown in Fig 10b1 (schb), b2 (sulpb), and b3 (bothb). Indeed, the attractor in Fig 10a1–10a3 could be morphed into the one shown in Fig 10b1–10b3 by twisting the upper half of the attractor with respect to the lower half.

As mentioned in the lag time section, the autocorrelation-based lag time (see Table 3) was an order of magnitude larger than the AMI-based lag time (see Table 4). This raised the question whether such low lag times could properly unfold the three dimensional attractor as we observed for much larger autocorrelation-based lag times in Fig 10. We show below a side-by-side example of an attractor for the same condition and the same trial unfolded both with the autocorrelation-based and AMI-based lag time in Fig 11. For the example shown in Fig 11, the condition, the mouse, and the trial were randomly selected from among the six conditions, 24 mice, and 100 trials for each mouse.

Fréchet distance in neural activity space

Finally, we performed a third statistical test of similarities among trials based on the reconstructed tree-dimensional attractors. For this purpose, we computed the Fréchet distance [133, 135] between three-dimensional phase space trajectories, such as those shown in Fig 10. Intuitively, Fréchet distance relates to the “man walking a dog on a leash” paradigm, i.e. find the shortest leash needed such that the man walks along one curve and the dog walks along the other [134]. Fréchet distance is a topological measure of the similarity between two curves [134], i.e. it does not consider any temporal constrains regarding the location of figurative points on the two curves. We computed Fréchet distance using the algorithm in [136]. In Fig 12, each white-bordered square corresponds to a combination of 100 baseline trials × 100 trials with DA receptor antagonists. For example, the top left square shows the color-coded Fréchet distance between the 100 baseline trials bc8 and the 100 trials after DA receptor antagonist sulpb8.

Fréchet distance is a topological measure of similarity between two curves, i.e. once the lag time was used for the purpose of unfolding the attractor the structure is “frozen” in phase space. This means that Fréchet distance only computes the topological distance between any two points of the two trajectories regardless of the speed at which the figurative points travel along the path, i.e. disregards dynamical aspects.

The fact that two different conditions, such as bc8 and sulpb8 have similar-looking three-dimensional attractors at an average minimum distance of 0.03 ± 0.01 arbitrary units (indicated by blue color in Fig 12) means that the figurative points of the two conditions visit similar phase space regions of neural activity. At the same time, for the same two conditions the weighted average lag times is 0.244 s both for the baseline (Fig 6b2) and sulpb8. This means that the two trajectories are similar both geometrically and that they are traveled at the same speed. Other conditions, such as bc9 and sulpb9, have quite large Fréchet distance of 0.15 ± 0.02 between reconstructed attractors despite the fact that their lag time distributions have the same weighted average lag times. The implication is that although the time constants for the phase space dynamics are identical for mice # 8 and 9 both for the baseline and DA receptor antagonist, the differences in topological distances suggest that either the initial conditions for the two mice were different or the sensitivity to light stimuli for the two mice is systematically different.

Towards a quantitative model of the experimental data

The reconstructed attractors (see Figs 10 and 11) suggest closed phase space, i.e. periodic, trajectories. While in principle any signal could be represented as an infinite series of harmonics, we only decomposed the signals using 3, 5, 10, and 100 of the harmonics with the highest amplitudes (see Fig 13). As the inset in Fig 13a1 shows, including more frequencies in the spectral analysis leads to a better approximation of the original LFP. The distribution of residuals, i.e. the difference between the LFP and the synthesized signal, shown in Fig 13a2 suggests that the spread of the distribution decreases as the number of harmonics increases. When 100 frequencies with the highest amplitude are considered (Fig 13b1), then the synthesized signal can reproduce minute details of the LFP and the corresponding residual has a Gaussian distribution with a very narrow range (see Fig 13b2).

We tested the LFP residuals obtained by subtracting the synthesized signal according to Eq 2 from the original LFP data (see Fig 13a2 and 13b2). In all signals synthesized with more than 10 frequencies the residuals were identified with Matlab system identification toolbox as ARMA(0,0) processes, i.e white noise. Additionally, we ran the augmented Dickey-Fuller and Kwiatkowski, Phillips, Schmidt, and Shin (KPSS) tests. Both tests showed that the residual time series shown in Fig 13 were stationary. When we used less than 10 frequencies of the Fourier spectrum to synthesized the signal, the residual signal failed one or both stationarity tests. This shows that with a too low number of harmonics the residual signal still contains trends that were not removed from the original signal when subtracting the synthesized one.

The summary of spectral analysis over all trials of sulpa for mouse #18 is shown in Fig 14a1 (first three strongest frequencies in Fourier spectrum), a2 (first five frequencies), and a3 (first ten frequencies). As Fig 14a1–14a3 suggest, for each frequency in the Fourier spectrum (horizontal axis) the corresponding amplitude has a consistent distribution (see the insets) regardless the number of selected frequencies. Fig 14b1–14b3 summarize the distribution of Fourier amplitudes over all trials of sulpa for mouse #18 for the first four frequencies.

We present in the Supporting Information the extensive S1 Table that contains all the frequencies and Fourier amplitudes necessary for the synthesis of the signal with three, five, and 10 highest amplitudes. We did not include in the table the model with 100 highest amplitudes since it only adds extra coefficients with little extra information compared to the already presented reduced models. First, we noticed that different trials may have different spectral decompositions. For example, when decomposed in three Fourier components, trial #1 requires frequencies k = 2, 4, 6 whereas trial # 15 requires frequencies k = 2, 4, 8. As a result, for each of the three decompositions in Fourier components, we also show what fraction of the trials required a particular frequency. Based on the number of data points and the sampling duration, the actual frequency is f_k = k/3 (Hz). Since in Fig 14 we only showed the distribution of amplitudes for the first four strongest frequencies, we show below a very small excerpt from S1 Table.

We notice from Table 5 that the first four frequencies have stable amplitudes over different Fourier representations of the signal, which is an indication that they capture a significant amount of information regarding the LFPs. Even in the case of a flexible decomposition of LFPs with 10 frequencies, the four strongest harmonics listed cover over 44% of the trials and 61% of the trials when constraining the Fourier decomposition to five frequencies (see S1 Table for detailed information on Fourier coefficients and their corresponding weights in data modeling).

To conclude, a spectral model of the data with less than ten harmonics is accurate enough. Given the large number of trials that are captured correctly even by only four harmonics (see Table 5), a much lower number of Fourier coefficients (maybe five) may suffice.

Discussion

What did we learn from analyzing the neural activity under DA receptor antagonists? First, we found that the delay embedding method works (see [102] for a discussion of control data and [103] for cocaine data). Second, we found that the embedding dimension is still three, the same as in the previous studies [102, 103]. This means that the mathematical model that could describe the control, the cocaine, and DA receptor antagonists only require three independent variables.

We also found that there are mice for which the overlap between baseline and DA receptor antagonists is maximum, i.e. 50% and, therefore, the dendrogram with Euclidian distance measure cannot distinguish between baseline and DA receptor antagonists. Such examples are mice # 11, 12, 20, 21, and 23. We also caution against exclusively and solely relying on dendrogram for the purpose of trials’ classification and statistics. For example, the dendrogram is sensitive to the type of distance selected for classifier.

Weighted average of lag times

The second, more qualitative classifier was based on the distribution of lag times. The delay times for the phase space reconstructions significantly differ among animals (see Fig 5). Although the optimal delay time even differs from trial to trial, it is clear that the distributions of delay times are also significantly different between baseline (Fig 6) and DA receptor antagonists for most mice (see Fig 5).

For some mice, the LFPs of baseline and DA receptor antagonists were not distinguishable using dendrogram method, whereas the lag time statistics allowed a clear separation of the two conditions. For example, for mouse # 12 the weighted average of lag times for baseline was 0.297 s which shifted to a lower lag times of 0.26 s when SCH 23390 was applied before cocaine (see Fig 7b1).

On the other hand, the weighted average of lag times for other mice did not change between baseline and DA treatments. For example, for the mouse # 11 the weighted average of lag times was 0.244 s (Fig 5b2), for mice # 20 and 21 it was 0.255 s (see Fig 6b3), and for mouse # 23 it was 0.178 s (see Fig 6a3). At the same time, only comparing DA receptor antagonists against baseline does not cover the entire complexity of the experiment we carried out. This is because the common denominator of all DA receptor antagonist experiments in this study is cocaine: we either injected cocaine before or after the D1 receptor antagonist SCH 23390, or the D2 receptor antagonist sulpiride, or both SCH 23390 & sulpiride. Therefore, we also checked DA receptor antagonists’ effect on lag time distributions when compared against pure cocaine experiments [103]. We are again interested in mice # 11, 12, 20, 21, and 23 since the dendrogram could not distinguish between baseline and DA receptor antagonist trials in these cases. As it turned out, in all of the above cases there were significant differences between the lag times of cocaine trials (see [103]) and DA receptor antagonists. For example, the weighted lag time for mouse # 11 is 0.285 s under cocaine and 0.244 s when sulpiride was applied before cocaine (Fig 5b2). The weighted lag time for mouse # 12 was 0.263 s under cocaine and 0.26 s when SCH 23390 was applied before cocaine (Fig 5b1). The weighted lag times for mice # 20 and 21 were 0.219 s for cocaine and 0.255 s for sulpiride and SCH 23390 applied before cocaine (see Fig 5a3). The weighted lag times for mouse # 23 was 0.242 s after cocaine and at 0.178 s for both sulpiride and SCH 23390 after cocaine (see Fig 5a3).

It is also important to note that when compared against baseline, the only trials in which the weighted average of the lag time increased were scha and sulpa. In all other conditions, i.e. schb, sulb, and both, the weighted average of the lag times remained the same. This indicates that from a dynamic point of view, it matters in which order is the SCH 23390 and sulpiride, respectively, injected with respect to cocaine. In trials where the DA receptor antagonists were injected after cocaine the lag time increased, which means a longer correlation time of the LFP recordings and, therefore, a different time scale of the model equations compared to the baseline. At the same time, when cocaine was applied before DA receptor antagonists we did not observe changes in the weighted average of the lag times compared to baseline. To fully understand these cases we need to compare them against the effect of pure cocaine (no DA receptor antagonists). We noticed that in sulpb, such as sulpb9 and sulpb11, the cocaine alone decreased the weighted average of the lag times compared to baseline whereas sulpiride applied before cocaine shifted the lag times to longer durations and precisely back to the baseline, practically canceling the effect of cocaine. Similar responses were observed in schb13 and schb15 where cocaine alone decreased the weighted average of the lag times compared to baseline and then SCH 23390 applied before cocaine practically canceled the effect of cocaine leaving the lag time at baseline values. Finally, in all cases when both SCH 23390 & sulpiride were applied before cocaine, i.e. bothb20, bothb21, and bothb22, the cocaine alone systematically moved the weighted average of the lag times to shorter durations compared to baseline, whereas the two DA receptor antagonists together practically canceled out the effect of cocaine and left the lag times at baseline values.

It is also important to address the issue of some trials that do not fit the above picture, such as sulp8 and sulp10 for which cocaine alone shifted the weighted average of the lag times to longer durations. Even in these two cases, the effect of sulpiride was consistent with the other sulpb cases, i.e. sulpiride applied before cocaine shifted the lag times back to baseline. The mismatch does not concern the effect of DA receptor antagonists, which consistently cancel out the effect of cocaine when applied before it, but rather the unusual increase in the weighted average of the lag times under cocaine. The cause could be our choice of the weighted average of the lag times as a measure for comparing two distributions, i.e. baseline and cocaine. In contrast, in [103] we first fitted the distributions with a log-normal function and then compared the location of the peaks of the distributions. In particular, sulpb8 and sulpb10 have a longer tail compared to sulpb9 and sulpb11 that shifts the weighted average of lag times towards longer durations while the peak of the log-normal distribution remains unchanged (not shown). Therefore, the mismatch could be determined by our choice of the number we put on a distribution of lag times. Furthermore, for the particular mice # 8 and 10, a log-normal fit is not quite optimal as the lag time distributions seem to have bimodal shapes.

Kolmogorov-Smimov (KS) similarity of lag time distributions

The KS test that compared the actual distributions of baseline lag times against DA receptor antagonists (see Fig 9b1) has a significantly large number of diagonal elements, i.e. cases where KS test placed the baseline and DA receptor antagonists in the same empirical distribution. Among others, we find again mice # 11, 12, 20, 21, and 23 in this category. Interestingly, DA receptor antagonists removed mice # 16, 17, 18, and 24 from the cross-list and indicated that their lag time distributions changed compared to baseline (see blue solid squares in Fig 9b1). This confirms the findings based on the weighted average of the lag times and supports the conclusion that DA receptor antagonists applied after cocaine, i.e. scha16, scha24, sulpa17, sulpa18, increased the lag time/correlation time of LFP recordings. The two missing mice for DA receptor antagonists applied after cocaine, i.e. sulpa19 and sulpa23 had a large p-value, i.e. close to the borderline of the statistical test.

As with any statistical test, we must also be aware of their weaknesses. For example, the KS test may not give reliable estimates for similar but different cumulative distribution functions since this statistical test is most sensitive near the median of the samples and less sensitive near the tails [146]. As we notice from Fig 5, all our lag time distributions are skewed towards one (very low) extreme of the spectrum. Therefore, we expect that KS test may fail on some of these distributions.

Three-dimensional reconstructed attractors

We reconstructed three-dimensional attractors based on LFPs from mPFC of ChR2 expressing PV+ interneurons (see Fig 10). We found topologically similar phase space attractors that could be morphed into each other after an appropriate change in delay time (see [102, 103]) for details). The characteristic “8” -shaped attractor (see Fig 10) was also found in the control [102] and cocaine cases [103].

All studies on control data [102], cocaine [103], and this DA receptor antagonists study suggest that the local network could be described by a model with only three variables. Furthermore, all control, cocaine, and DA receptor antagonists trials are classified as similar and placed in the same cluster by the dendrogram method in at least 20% of the cases, regardless of the number of clusters of the dendrogram (see Fig 4c). This overlap suggests that there may be a common, invariant, part of the mathematical model that describes all trials. At the same time, there are some significant difference in the delay times between control and DA receptor antagonists (see Fig 5) that could lead to different phase space dynamics. The differences in the delay times between baseline and DA receptor antagonists, whenever significant, could be interpreted as two different time scales for the two experiments.

We also found that the attractors correctly unfold for a wide range of lag times. Both the attractors unfolded with the autocorrelation-based lag time (see Fig 11a) and with the AMI-based lag time (see Fig 11b) are correctly unfolded. The only difference is that the attractors for unfolded using the very short lag times obtained with the AMI-based method are very narrow along a diagonal.

Fréchet distance between pairs of attractors

Fechet distance plots shown in Fig 12 provide an intuitive and global understanding of potentially similar time scales between different trials as represented by blue colors. Although the 3D attractors use the lag time to unfold the phase space trajectories (see Fig 10), they lack any temporal dimension. Therefore, Fréchet distance is a measure of only topological similarities between phase space trajectories. In other words, Fréchet distance can tell how close the attractor visits a certain region of the phase space, but not how often or how long the system remains in that region of the phase space. As a result, Fréchet distance offers a complementary view of reconstructed dynamics when compared against lag times, one in which we are interested more in the geometry rather than the dynamics of the phase space. One very good example of the complementarity is bc8 attractor that is very similar topologically to both sulpb8 and, surprisingly, sulpa19 (see Fig 12). From a dynamic point of view, the KS test of lag/correlation times distributions for bc8 and sulpb8 shows that they were identical (Fig 9b1), but bc8 and sulpa19 do not have similar lag time distributions. This suggests that different statistical measures capture different aspects of neural activity, some capture the time scale of the dynamics (lag/correlation time) and others capture the topology of the phase space trajectories and the type of accessible states.

Fourier analysis

The closed phase space trajectories shown in Figs 10 and 11 suggest that spectral (Fourier) analysis may capture the most significant part of the LFP dynamics. We found that Fourier decomposition with more than ten harmonics leaves a white noise residual. Depending on the required accuracy of the model, it is even possible to capture the main trends with less than five Fourier coefficients.

Conclusions

The activity of the medial prefrontal cortex (mPFC) in 17 mice across six difference combinations of two DA receptor antagonists were investigated. The D1 receptor antagonist SCH 23390 and D2 receptor antagonist sulpiride were used in these experiments either applied before, as in schb, sulpb, and bothb trials, or after, as in scha, sulpa, and botha trials. The transient phase resetting induced by a light stimulus was removed using the pair correlations between recorded local field potentials (LFPs) [102, 103]. While the improvement in the correlation coefficient is qualitatively apparent from Fig 2, the compact numerical summary in Table 3 shows that the circular shifting improvement exceeds one order of magnitude.

The phase-corrected trials were embedded in a three-dimensional phase space using a delay-embedding method. The main novel results of this study are as follows:

• We reconstructed the DA receptor antagonist attractor and found that it is three-dimensional.

• In all trials SCH 23390 after cocaine increases the lag time by about 55% compared to pure cocaine case and when applied before cocaine it decreases the weighted lag time by about 1%. The same trend is preserved when comparing SCH 23390 against baseline, albeit at different percentage changes (see Table 3).

• Sulpiride after cocaine increases the lag time by about 26% compared to pure cocaine trials, whereas when applied before cocaine it decreases the weighted lag time by about 14%. Sulpiride after cocaine preserves the trend when compared against baseline, although when applied before cocaine produces no change compared to baseline (see Table 3).

• When both SCH 23390 & sulpiride are applied they seem to reverse the trend observed in the single DA receptor antagonist cases above. For example, when both SCH 23390 & sulpiride were applied after cocaine they seem to decrease the weighted average of lag times by 27%, although the observation is only based on one available mouse. When both SCH 23390 & sulpiride are applied before cocaine its weighted average of lag times increases by 16%. There is no change in the weighted average of the lag times between baseline and both SCH 23390 & sulpiride receptor antagonists case.

• A spectral model for the LFP recordings with as low as five Fourier coefficients can capture a significant portion of data trends.

To summarize, both SCH 23390 and sulpiride, either applied separately or together, cancel out the effect of cocaine when applied before it (see Table 3). At the same time, both SCH 23390 and sulpiride, either applied separately or together, consistently increase the lag/correlation time of LFP recordings when applied after cocaine.

Supporting information

S1 Fig [avi]
Animated three-dimensional reconstructed attractor using autocorrelation-based lag times.

S2 Fig [avi]
Animated three-dimensional reconstructed attractor using AMI-based lag times.

S1 Table [hz]
Fourier coefficients for signal synthesis.