Cocaine diminishes functional network robustness and destabilizes the energy landscape of neuronal activity in the medial prefrontal cortex

Abstract We present analysis of neuronal activity recordings from a subset of neurons in the medial prefrontal cortex of rats before and after the administration of cocaine. Using an underlying modern Hopfield model as a description for the neuronal network, combined with a machine learning approach, we compute the underlying functional connectivity of the neuronal network. We find that the functional connectivity changes after the administration of cocaine with both functional-excitatory and functional-inhibitory neurons being affected. Using conventional network analysis, we find that the diameter of the graph, or the shortest length between the two most distant nodes, increases with cocaine, suggesting that the neuronal network is less robust. We also find that the betweenness centrality scores for several of the functional-excitatory and functional-inhibitory neurons decrease significantly, while other scores remain essentially unchanged, to also suggest that the neuronal network is less robust. Finally, we study the distribution of neuronal activity and relate it to energy to find that cocaine drives the neuronal network towards destabilization in the energy landscape of neuronal activation. While this destabilization is presumably temporary given one administration of cocaine, perhaps this initial destabilization indicates a transition towards a new stable state with repeated cocaine administration. However, such analyses are useful more generally to understand how neuronal networks respond to perturbations.


Introduction
The intricacies of a brain's neuronal network drive its functionality.The characterization of such a network spans multiple scales ranging from the scale of genes to the scale of the whole brain.For instance, specific genes have been identified as playing a role in brain size and shape and, more recently, in long-range connectivity (1)(2)(3)(4)(5)(6).At the brain scale, on other hand, there exists the nearly universal finding of a small-world network architecture in which most nodes can be reached by traversing some number of edges (7,8).Others focus on particular neuronal circuits at the mesoscale to decode a specific functionality.For instance, it was recently revealed how a neuronal circuit involved in fruit fly navigation empowers a fly to perform vector addition (9).How such findings at the different scales can be integrated into one, more complete picture to understand the structure-function relationship in the brain is a current avenue of investigation (10,11).
At the mesoscale, in brains with smaller numbers of neurons, such as the Caenorhabditis elegans or the fruit fly, determining neuronal circuits is feasible.In brains with larger numbers of neurons, determining such neuronal circuits is not as feasible.And yet, with the advent of new imaging techniques in freely moving animals, such as rodents, to record individual neurons deep inside the brain (12), one is compelled to assess the functional contribution of a collection of individual neurons forming a subnetwork.Here, we will study such a subnetwork and characterize its response to a particular perturbation.Our perturbation of choice is the impact of cocaine on a neuronal subnetwork of the medial prefrontal cortex.
Cocaine primarily impacts the brain's limbic system by targeting dopaminergic pathways and increasing dopamine release (13,14).The limbic system is a network of multiple brain regions including the prefrontal cortex and the nucleus accumbens, which are brain regions closely associated with motivated behavior and substance use disorders.The prefrontal cortex acts in a top-down manner to inhibit the nucleus accumbens and suppress motivated behavior (15).Therefore, stable function of the prefrontal cortex is critical to suppressing the abberant motivated behavior associated with substance dependence.
And so while the interconnectedness of the brain's limbic, or emotional, system is key to understanding the effects of cocaine on the brain, here we focus on the medial prefrontal cortex.At the scale of the prefrontal cortex, there is a reduction in the prefrontal cortex volume in cocaine-dependent men that may help explain deficiencies in prefrontal cortex functionality (16).Moreover, electrophysiology has shown that the administration of cocaine decreases neuron activity in the medial prefrontal cortex (17).Additionally, from a neuronal network point of view, the likelihood of developing dependence is associated with significant hypoconnectivity in orbitofrontal and ventromedial prefrontal cortical-striatal circuits-pathways critically implicated in goaldirected decision-making (18).At smaller scales, there exists some evidence of neuronal remodeling at the single neuron level in the medial prefrontal cortex (19)(20)(21).
The impact of cocaine on the medial prefrontal cortex, in particular, has relied heavily on functional magnetic resonance imaging to understand the function of the region as a whole (16,22,23) or in vitro or ex vivo analyses of gene or electrophysiological profiles to understand single neuron alterations (19,20,24).Our approach, on the other hand, utilizes a combination of established and novel statistical mechanics methods to identify cocaine-induced alterations in the prefrontal cortex network at the mesoscale.Furthermore, we provide a new analytical tool for in vivo calcium imaging data from freely behaving rodents as well as show, for the first time to our knowledge, destabilizing effects within a statistical mechanics setting, of single-cocaine administration on the prefrontal cortex.
We will focus in on the mesoscale by recording from a subset of neurons in the medial prefrontal cortex of a living rodent before and after cocaine administration.In particular, we will model the neuronal subnetwork as a linearized version of a Hopfield model (25,26) and use a machine learning approach to obtain its functional connectivity.Determining the functional connectivity map allows one to make powerful predictions about the dynamics of the network in terms of firing patterns in response to external inputs, for example.In addition, it also allows one to perform network analysis (27,28) with the additional information of functional-excitatory and functional-inhibitory weights that are typically not known when analyzing a structural connectome from a network point of view (29).Many efforts are indeed underway to correlate the functional connectome with the structural connectome (30).On the other hand, even in the absence of such information, network analysis can identify subgraphs, which may correspond to a neuronal circuit with a particular function (31).
At the mesoscale, one can also toggle between dynamical approaches and statistical mechanics approaches, in which time-averaged information may be useful.Dynamical approaches with large groups of neurons are less feasible, at least in so far as obtaining the functional connectivity of the group.As for prior statistical mechanics approaches, much has been done to give rich, new insights into brain functionality such as Refs.(32)(33)(34).Moreover, statistical field theory approaches to neuronal networks also exist (35)(36)(37).We will use a mean field theory approach in which spatial degrees of freedom are averaged out to obtain new insights into a Landau-Ginsburg-like energy landscape for the effects of cocaine on the brain.
The outline of the manuscript is as follows.First, we detail the experiments.Then, we briefly outline both the theoretical framework of the Hopfield model, applied to a neuronal network, as well as a statistical mechanics approach.After presenting the results of the analysis from both approaches, we conclude with a discussion of the impact of our results.

Animals
Male Sprague Dawley rats (n = 3) (Harlan Inc, Houston, TX, USA) aged 10 weeks were initially pair housed and maintained on a diet of standard rat chow (Tekland Mouse/Rat Diet 7912, Harlan Laboratories, Inc., Indianapolis, IN, USA) ad libitum in their home cage.The colony was maintained on a standard 12-h light cycle (6:00-18:00), 71 • F with relative humidity of 30-50%.All animal use was carried out in accordance with the Guide for the Care and Use of Laboratory Animals and with the approval of the University of Texas Medical Branch Institutional Animal Care and Use Committee.

Virus injection
Following a 1-week habituation period, rats were anesthetized using 1-5% isoflurane and placed in a stereotaxic surgical apparatus (Kopf, Tujunga, CA, USA).Each animal received a unilateral injection of 1.0 μL AAV1.Camk2a.GCaMP6m.WPRE.SV40 (Inscopix, Mountain View, CA, USA) at a rate of 0.2 μL per minute in the right hemisphere medial prefrontal cortex (AP + 3.0 mm, ML + 0.5 mm, DV −4.2 mm).The syringe remained in place for 5 min following completion of injection before being removed and the wound was closed.Rats were then placed in a clean cage on a heating pad and followed postsurgically to ensure recovery.Animals were habituated to daily handling following surgery.

Gradient-index lens implantation
Four to five weeks after virus injection, rats were anesthetized using a 1-5% isoflurane.Three stainless steel surgical screws (one per skull plate) were inserted to provide stability for the lens.After screws were placed, rats underwent stereotaxic implantation of the ProView integrated gradient-index (GRIN) lens (0.6 mm × 7.3 mm, Inscopix) (AP + 3.0, ML +0.5 mm, DV −4.00 mm) in the medial prefrontal cortex.Excess space between the craniotomy and the base of the lens was filled with adhesive cement (Metabond Quick!Adhesive Cement System, C&B), and then the lens was secured to the skull with dental cement (ACRAWELD Repair Resin, Henry Schein, Melville, NY, USA).A baseplate cover (Inscopix) was installed to protect each lens until imaging.Wound clips were applied, and the animals were placed in a clean cage on a heating pad and monitored postsurgically to ensure recovery.Following GRIN lens implantation, all animals were single-housed to protect the implant and were handled daily.Postmortem histology was done to verify targeting.

Cocaine administration and calcium data acquisition
Once animals had recovered from lens implantation surgery (>3 weeks), they were habituated to light restraint and attachment of the miniscope (Inscopix).The miniscope was mounted and animals were placed in the testing chamber for recording at least three times prior to testing.The testing chamber was a 40 cm × 40 cm clear, Plexiglass container with cornbob bedding.During these recording sessions, field of view and imaging settings were optimized for each animal.All recordings were taken at 10 frames per second.On testing day, animals were given 5 min following miniscope attachment to acclimate to the chamber before recording began.Imaging for each animal was optimized, with gains 2.0-3.5 and LED power 0.5-1.3.Focal depths for each animal were chosen to capture the maximum number of clearly defined neurons in a single plane of recording.Thirty minutes of baseline calcium activity was recorded, followed by a 5-min pause during which animals received a single intraperitoneal injection of 15 mg/kg cocaine (National Institute of Drug Abuse, Research Triangle Park, NC, USA) in sterile 0.9% NaCl before being returned to the chamber.Thirty minutes of postadministration calcium data was recorded.Then the session was ended, and the animal was returned to its home cage.

Calcium activity trace extraction
Initial calcium imaging data processing was done using the Inscopix Data Processing Software (Inscopix).Video processing steps were the same for every animal.Videos were spatially downsampled by a factor of four.Motion correction utilizing the "mean image" as a reference frame was done and CNMFe was used to identify neurons (38).Manual validation of CNMFe identified neurons was performed to ensure that all identified neurons used for analysis contained a clear baseline for the 30 min and welldefined, positive peaks of fluorescence.See Fig. 1.The calcium traces are, therefore, ΔF/F, where F denotes relative fluorescence intensity with respect to the baseline.

Mathematical framework
In this section, we lay down the mathematical basis for the analysis we will use to determine the effects of cocaine on the medial prefrontal cortical network.We start from the modern Hopfield model of a neuronal network to write the governing differential equations of neurons' firing states, and to construct the total energy of the system (26).We review how the classic Hopfield model can then be extracted from the modern version (25) as well as address its biological interpretation using a circuit analogy (39) for a simpler version.Next, we construct the statistical field theory of a neuronal network starting from the Hopfield model energy.Within the context of both approaches, we discuss how the firing behavior of the network can be influenced by external stimuli.

Hopfield model of a neuronal network
The modern Hopfield Network (26) provides a mathematical description of the activity of neurons in the brain, which account for functional features, such as memory.In this model, neurons that are categorized as the feature neurons are referred to by Latin indices and hidden neurons are referred to by Greek indices.The feature neurons are the ones observed in an experiment and the hidden neurons are the rest of neurons that actually exist in the brain's region of interest.The two categories of neurons are linked by a matrix ξ iμ that represents the synaptic connectivity.The currents of the two types of neurons are denoted by V i and h μ .The output currents of the feature and hidden neurons are denoted by g i and f μ and are defined as where L h and L V are the Lagrangians that define the model.Therefore, f μ and g i are the conjugate variable to the currents of hidden and feature/visible neurons, respectively.The full theory is described by two differential equations given by where V (ext.)i and h (ext.) μ represent external currents that do not come from the rest of neurons in the system, τ is a time constant, and the sums run over all the neurons.The Hamiltonian of the neuronal network that drives its time evolution is given by ( Assuming that (i) the hidden neurons' activities are rather fast and (ii) they do not receive external currents, we can set τ h = h (ext.)μ = 0 and solve Eq. 2 to remove the current of the hidden This solution to the hidden current can be substituted into f μ in the first line of Eq. 2 such that the differential equation of the current of the feature neurons is self-contained.At this point, we need to specify the form of the Lagrangian, and, therefore, the Hamiltonian, to remove f μ in terms of V i and ξ iμ .We choose a very simple form, namely, ( This form leads to a simpler linear version of the classic Hopfield model (25).As we will observe in the next section, this simpler linear version, as opposed to piecewise linear, readily approximates the calcium intensity signal, which has a longer time constant than action potentials.Given the longer time constant, combined with experimental limitations, observing thresholding in individual neurons is difficult and so a simpler integrate-and-fire neuron model for f μ is more relevant here (40).
As has been done with the classic Hopfield model ( 25), as well as other neuronal network models (40), it is helpful to construct a circuit diagram of the neuronal network.We focus on one neuron in the network.More precisely, (i) a given neuron is replaced by a capacitor C with a resistor R in parallel, (ii) the remainder of neurons are replaced by a battery, and (iii) there exists a switch S that alternatively connects and disconnects the battery to the given neuron.Should the sum of functional-inhibitory and functionalexcitatory input signals for a given neuron is above a threshold, the battery is connected to the given neuron to lead to a charging phase followed by discharging phase in which the neuron does not fire regardless of the input signal from the rest of neurons.When the given neuron is disconnected from the battery S, it forms a self-loop.The model is presented in Fig. 2 on the left.When the battery is connected, the electric current across the membrane of the neuron satisfies the following differential equation where τ = RC with R being the resistance and C being the capacitance, the subscript i refers to the neuron that is modeled as an RC circuit, and i is the current of the battery.Index j refers to all the rest of neurons while T ij represents the functionalexcitatory or functional-inhibitory signal that neuron j sends to the neuron i.Using Eqs. 2 and 5, T ij has the following form It is this matrix, the functional connectivity matrix, that we will extract from the experimental data to forecast the neuronal network's future, assuming there is no plasticity over some time scale.
Right after neuron i is fully charged, the discharge state begins when the switch S disconnects the battery and makes a self-loop where the neuron satisfies the following differential equation Solving the two differential equations above, we readily find the current of neuron i in a full cycle whose solution is graphically depicted in right figure in Fig. 2 and given by where T refers to the time that the switch disconnects the battery.

Statistical properties of the neuronal network
It is also interesting to look beyond the microscopic details of the system, in terms of a functional connectivity matrix, to quantify its macroscopic properties as well in a statistical sense.In other words, the emergent properties of the neuronal networks can also be understood through studying the macroscopic behavior of the system using techniques from statistical mechanics as has been done previously (32)(33)(34)(35)(36).
We now consider fluctuations of the neuronal subnetwork, due to the surrounding neurons and other smaller scale influences.As this living system is not in equilibrium, we invoke a nonequilibrium version of equilibrium statistical mechanics where T eff is some effective temperature (41)(42)(43).Moreover, the unnormalized probability of occurrence of a state with energy E, after integrating out the hidden neurons in Eq. 3, is equal to e −E , with the k B T eff = 1, as it maximizes the entropy of the system (44)(45)(46)(47).Therefore, the partition function, i.e. the weighted sum over all the energy states reads where N ′ is the normalization factor.
We will assume that the energy fluctuates around some local minimum E 0 such that where  V = (V 1 , V 2 , . . . ) encodes the fluctuations of the system, the external source to the system of neurons is represented by  J, the Greens' function G is an unknown to be derived from data, and the rest of the terms that do not depend on V are denoted by E no−V .It is assumed that there exists at least one stable, local minimum.Of course, there may be multiple local minima that could complicate our analysis.The partition function can then be rewritten in the following form where the new normalization factor N , the path integration  DV, and the energy In principle, we should be able to use the partition function in Eq.
11, together with  V, and study the collective behavior in this space.
In practice, due to the high dimensionality of  V and the typically small number of collected data points, we choose to define a different parameterization with lower number of features.In other words, let us use a mean field representation for the average of all the neurons firing, or neuronal activity, at a given time, i.e.
where V 0 is some reference voltage.Therefore, we reorganize the sum in the partition function in Eq. 11 and define a new energy where  DV| m refers to a summation on all the possibilities in  V space that are constrained to a fixed m as defined in Eq. 13.Also,  Dm refers to the summation over all the values of m.Finally, given the above expansion about a minimum, we consider a similar expansion in m, or where a i are constants to be determined by data.Since these parameters determine the collective behavior of the neuronal network in m space, or the mean field neuronal activity space, we would like to learn see their variations by changing the experimental condition.We will address this possibility in the Results section.Interestingly, there has been a study quantifying the mean activity of the system as a function of time using the Wilson-Cowan equations (48,49) for neuronal activity as a starting point and then performing a similar expansion in the neuronal activity (37).Here, we investigate the distribution of the mean activity over time to assess its form.Moreover, as we will see, the data are well characterized by the expansion.

Results
We now use the data collected from the medial prefrontal cortex before and after cocaine injection as described in Experiments section to find the free parameters of the mathematical model in Mathematical framework section.First, we extract the synaptic connectivity of the neurons for each of the two experimental conditions.Next, we use the statistical field theory approach to predict the firing behaviors of the networks of each experimental conditions when they receive planned external stimulus.

Inferring the functional connectivity of the neuronal network from neuronal activity
Here, we derive the T ij connectivity matrices, and V (ext.)i in Eq. 6 by using a machine learning based model for each of the two datasets of Experiments section.From now on, we refer to the two experimental conditions as "before the injection" and "after the injection."It should be noted that we use a CamKII promoter for GCaMP expression.While this promoter has been thought to be specific to excitatory glutamatergic neurons, multiple lines of study show this specificity is not present, particularly in cortical circuits (50,51).As such, we make no assumptions about the biological identity of individual neurons and use functionalexcitatory and functional-inhibitory to describe the directionality of the inferred correlations between neurons.Therefore, the same neuron can participate in both functional-excitatory and functional-inhibitory relationships within the overall network.Stated another way, a neuron's calcium fluorescence may be positively or negatively correlated with the fluorescent intensity of other neurons in the network, which we describe as functional-excitatory or functional-inhibitory, respectively.
To remove background noise from each neuronal trace, each neuron, measurements of relative calcium intensity with a value smaller than a threshold are set to zero.We determine this threshold using the following method.We perform a sensitivity analysis by varying the relative calcium intensity threshold parameter.The procedure involves systematically changing this threshold parameter within the range of (0-20) and observing how these changes affect the coefficients and intercepts of the time-series linear regression model to be described below.For each threshold value, we select a random 70% subset of data, fit the time-series linear regression model to this subset, and store the model's coefficients and intercepts.This process is repeated for each threshold value, generating a series of models.After all models are generated, we assess the stability of these models by comparing the coefficients and intercepts across successive threshold values.This comparison helps in understanding how sensitive the model parameters are to the changes in the threshold.The aim of this analysis is to identify a threshold that ensures model robustness, meaning the model parameters do not significantly change with small variations in this threshold.Fig. S2 shows the relative calcium ion intensity change of the first neuron in the set over the timing of the experiment.Fig. S3 shows the recorded relative intensities of the entire set of neurons and the applied selections prior and post cocaine administration.

Borzou et al. | 5
To begin to analyze the functional connectivity of the neuronal network, we first time discretize Eq. 6 assuming that the RC time constant is of order the time increment δt such that In this equation, we set δt = 0.1 s to be consistent with the frequency at which neuron activity is recorded, the time scale at which we explore the functional connectivity of the neurons.This equation has the same form as the formulation of the supervised auto-regression time-series forecasting machine learning method, where the value of a feature such as V i can be predicted using its values in the past.To perform the learning analysis, we prepare the predictor dataset X in a matrix form such that each component holds a value of V i in a way that the rows refer to the time of the feature in an ascending order and rows refer to each of the neurons.We also prepare the response variables y for each neuron separately.The latter is again the same V i values for the corresponding neuron but at one time step later.Having the X, and y variables, we use the LinearRegression model of Scikit-Learn package (52) in python to perform the regression and use the returned coefficients and intercept as the T ij and V (ext.)i of that neuron.We repeat the process for all the neurons.For each of the experimental conditions, we use 75% of the dataset for training the model, i.e. estimating T ij and V (ext.)i , and the remaining 25% for testing the quality of the extracted parameters.For the training part, we label the predicted value of Eq. 16 as Vi (t) and define the error function to be where N is the number of neurons, V i (t) is the observed value, and d runs over the data points in the training dataset.Our estimation for T ij and V , from the hidden neurons, to the feature neurons is also estimated through our minimization.We find that these are negligible with respect to T ij matrices.The external inputs are shown in Fig. S4.We then use the resulting estimate of T ij to understand how cocaine administration may alter network function in the medial prefrontal cortex.We see a significant change between the connectivity matrices T ij of before and after cocaine administration.The two matrices are shown in Fig. 3, where each of the 51 neurons in the experiment are shown on both x and y axes.The functional-excitatory and functional-inhibitory neurons, i.e. the positive and negative components of T ij , for both of the experimental conditions are shown in Fig. 3.A few conclusions can be readily drawn from the figure.First, the number of the functional-excitatory signals, shown by red colors, is close to two times the number of the functional-inhibitory signals, shown by the blue colors, in both of the experimental conditions.This result is in surprising agreement with biological observations of the actual synaptic connectivity (53,54).Second, the connectivity matrices T ij are not symmetric, which is again in agreement with neuroscience expectations (55)(56)(57).Third, the interneuron signaling prior to the cocaine administration, when the neuronal network is in its "normal" state, is less significant than the signaling activities post cocaine administration.However, interestingly, we observe that the changes in decreasing the functional weights and increasing them are somewhat symmetric (see Fig. 4).This finding suggests that, at least at this scale of 50 neurons, the activity of the medial prefrontal cortex is modified beyond an overall decrease in activity, as reported in electrophysiology measurements on the scale of the brain region (16).Our analysis captures more subtle changes in functionality of the neuronal network.It should be mentioned that the functional connectivity matrix is symmetric with 50 dimensions (neurons).Hence, its unique components total 1,275.Additionally, the external voltage contribution introduces another 50 variables.Consequently, to resolve all unknowns, a minimum of 1,325 data points is required.This translates to an approximate need of at least 2.2 min of data, taken at 0.1 s intervals, for a comprehensive analysis, providing a rough estimate for the lower bound on the amount of data required.To obtain further insight into the change in functional connectivity before and after cocaine administration, we performed conventional network analysis using the python package NetworkX.In NetworkX, the diameter of a graph is defined as the maximum eccentricity among all pairs of nodes in the graph.The eccentricity of a node is the maximum distance from that node to any other node in the graph.For weighted graphs, the distance between nodes is defined as the sum of the weights along the shortest path between them.The diameter function in NetworkX takes this into account when working with weighted graphs.Moreover, the betweenness centrality for a node v is defined as the fraction of all shortest paths, between all possible node pairs, in the graph that pass through v.For a weighted graph, the shortest path between a given pair is the sum over the weights of the edges along that path.
We find that for both the functional-inhibitory neurons and the functional-excitatory neurons that the diameter of each respective graphs increases after the cocaine injection.To be specific, the diameter of the functional-inhibitory graph is 0.31 before cocaine and 0.62 after.For the functional-excitatory neurons, the increase is even more significant with 0.46 before and 1.38 afterwards.A lower graph diameter typically correlates with a more robust graph in that it is a more tightly connected graph.Therefore, cocaine has diminished the robustness of both the functional-excitatory and functional-inhibitory network functionality.While these reported graph diameter increases are for one animal, the trend is similar for the remaining two animals.To be specific, the graph diameter of the functional-inhibitory graph is 0.18 before cocaine administration and 0.75 after, while for the functional-excitatory graph, the graph diameter goes from 0.32 to 0.39 for the second animal.For the third animal, the graph diameter of the functionalinhibitory graph is 0.24 before cocaine administration and 0.96 after; for the functional-excitatory graph, we compute a graph diameter changing from 0.46 to 1.93.
For betweenness centrality, after determining the shortest path between any two vertices, the betweenness centrality of a vertex is the number of such shortest paths that include that particular vertex.See Table 1.For the functional-excitatory effects, neurons 28 and 5 take on more of a role in the network with both of their betweenness centrality scores increasing.Both neurons are enhanced in terms of their capability to activate other neurons.Whereas for the functional-inhibitory effects, neuron 1 loses its ability to inhibit other neurons as its betweenness centrality score decreases to essentially zero, while the capability for neuron 17 (in the bottom 5 rows) to inhibit other neurons is enhanced.

Stimulating the neuronal network
A natural question after learning the connectivity map of the neuronal network is how we can control the firing pattern.Can we send external current to the network to create a certain neuronal firing pattern that leads to a specific function in the brain?
In order to see how an input signal to the network at time t affects the firing pattern at time t + n δt, we start from Eq. 16 and recursively apply it to derive the following equation where T is the matrix form of the connectivity matrix T ij without its diagonal components, and  V = (V 0 , V 1 , . . .).We send the same time-varying external current to all of the neurons.The external current increases for a while and falls to zero at some point.The left panel of Fig. 5 depicts the time dependence of the external current.The response of the neurons to the external input is shown in the right panel of Fig. 5.An interesting observation is that while all of the neurons received the same external current, the peak of their currents is not the same.Also, the decaying rate of the neurons' currents, after the external input is turned off, is not the same.More interestingly, the decaying rate and the peak of the currents are somewhat independent.For example, the neuron with the tallest peak decays faster than some other neurons.All of these observations can be explained by referring to the properties of the connectivity matrix T ij .The current in the neurons is induced by (i) the external input,  which is the same for all the neurons and (ii) the rest of the neurons.Therefore, the difference in the current of the neurons is a function of the signals they receive from within the network.For example, after the external current is shut down, the neurons continue to send signals to each other.Depending on their connectivity, they can receive signals for longer time, and as a result, their current decays more slowly than the rest.A similar plot is shown in Fig. S6 where two neurons are excited for a very short time and the current of the three most excited other neurons and three least excited other neurons are plotted.

Stability analysis in the collective neuronal activity space
While the prior subsection addresses how can use the functional connectivity to predict the firing pattern as a consequence of additional stimuli, we use our notion of energy in the collective neuronal activity space, or m-space, to study the changes in the parameters of the energy of Eq. 15 that are potentially induced by cocaine.First, we binarize the currents V i to 0 and 1 depending on whether they are below the threshold defined above.Next, we use the binary current to compute m at each given time.The time variation of m over the entire experiment, both before and after cocaine administration, is plotted in Fig. S5.The distribution of the mean-binarized-current m for before and after the injection is plotted in Fig. 6(left).
Our goal is to use the data in Fig. 6(left) to learn the free parameters of the energy, i.e. the parameters of the minus logarithm of the probability function.The popular method of maximum likelihood may not easily converge in this case.The reason is that the normalizing parameter a 0 in Eq. 15 is not analytically known due to the higher order terms of m in the probability function.
Our approach to estimate the free parameters of the probability function is to use the kernel method (58) to estimate the distribution of m from data in a nonparametric way.We scale the data, separately for each of the two experimental conditions, by subtracting the mean of m from each m and dividing by the standard deviation of m Next, we use the trained estimator to predict minus the logarithm of the probability over a grid of m, with one million bins, that spans from −2.7 to 2.7 in the scaled space.The range is chosen such that the probability is negligible beyond it.Next, we regress the one million estimations of minus the logarithm of the probability on  4 j=0 c j m , i.e. the Taylor expansion of the energy in the scaled space.Finally, after learning c j from the regression, we convert back to the nonscaled m space to find the values of a j in Eq. 15.
The learned parameters are shown in Table 2.
We use the learned parameters a j and Eq. 15 to find minus the logarithm of the probability states of m in before and after the injection.The results, shown in Fig. 6(right), indicate that the energy of the neural network is at its minimum prior to cocaine administration.After cocaine administration, the system is shifted to a state whose energy does not have a minimum.Therefore, prior to cocaine administration, the state of the neuronal network is stable.Soon after cocaine administration, the neuronal network is in an unstable state and may presumably evolve to go back to the stable state over time as effects of the cocaine subside.Hence, the connectivity matrix of Fig. 3(right), and its representation in Fig. 7(bottom), are transient in nature.We observe a similar destabilization of the energy in a second animal in response to cocaine (see Fig. S7).In a third animal, we observe an enhancement of destabilization in the energy (see Fig. S8).
Another interesting characteristic of Fig. 6(right) is that the energy post cocaine administration drops toward positive m, i.e. the probability of that direction is higher.Therefore, the system will evolve toward positive m direction.On the other hand, in Fig. 6(left), the mean of the prior to cocaine administration is to the right of the mean of before the injection data.Putting the pieces together, after cocaine administration, its neuronal network will presumably evolve back towards prior to cocaine adminstration m space.Finally, the shift to the left in m after cocaine suggests an overall average decrease in activity as consequence of cocaine.

Discussion
We present analysis of neuronal activity recordings from a subset of neurons in the medial prefrontal cortex of rats prior to and postadministration of cocaine.Using an underlying modern Hopfield model as a description for the neuronal network, we use a machine learning approach to determine the underlying functional connectivity of the neuronal network.The functional connectivity specifies which derived neuronal connections are functionalexcitatory and which neuronal connections are functionalinhibitory as well as the strength of the functional-excitatory or functional-inhibitory interactions.We find that the functional connectivity changes after the administration of cocaine with both functional-excitatory and functional-inhibitory neurons being affected.With such quantification, we can make predictions about the neuronal network, at least over some time period, in response to external stimulation of individual neurons or even collections of them.
To explore the impact of the changes in functional connectivity in response to cocaine, there are a number of different analyses one can perform.We perform conventional network analysis in terms of unweighted topological measures of the graph.We find that the diameter of the graph increases with cocaine, suggesting that the neuronal network is less robust.We also find that the betweenness centrality scores for a few of both the functionalexcitatory and functional-inhibitory neurons decrease significantly, while other scores remain essentially unchanged.Since betweenness centrality (for neurons) is a measure of the network throughput, the smaller the betweenness centrality on average, the more robust the network given that the shortest paths are more evenly distributed throughout the network as opposed to relying on a few nodes.The increase in betweenness centrality suggests that, again, the network has become less robust.
And, yet, while these measurements further quantify changes in the neuronal network in terms of graph theoretic measures, the analysis does not take into account individual neuronal activity directly.There exist methods to extract the functionality of the neuronal network, such as the controllability of neuronal networks (61).Since controllability analysis is more detailed with its various assumptions, we present a statistical analysis directly rooted in the activity of the network.Specifically, we have studied the spatially averaged neuronal activity over time and studied its statistics, just in the manner one would do with a spin system in which the sum of the neural activity relates to an order parameter of the spin system.Note that we are focusing on the neuronal excitation of the system, or the activation, and not the inhibition, as functional-inhibitory signals received by individual neurons cannot be measured with GCaMP.Just as in physical systems, a Fig. 6.Energy functional of neuronal activity.a) Distributions of the neuronal activity m in before and after cocaine administration.The plot indicates the differentiation in the collective behavior of the two systems due to the variation in the experiments.In this plot, the data prior to cocaine administration is weighed such that the two histograms have the same under-the-curve area.b) Minus the logarithm of the probability function, as defined in Eq. 15, learned from data in the left panel.This plot indicates that the energy of the neuronal network prior to cocaine administration is stable where the energy has a well-defined minimum.However, following cocaine administration, the energy of the neural network is pushed to an unstable state where the probability function does not have a minimum.nonequilibrium version of a Landau-Ginsburg-like theory.The shape of energy functional in m-space, as we denote it, tells one something rather important about the state space of the system.Our m-space analysis demonstrates that prior to cocaine administration, the neuronal network is stable in an energetic sense, while post cocaine administration, the neuronal network is unstable energetically.This discovery demonstrates that the administration of cocaine has strongly influenced the network to position it to be heavily influenced by other factors to be driven to a completely new stable state.In the language of substance abuse, the system is being readily driven to a new stable state that may underlie future dependence.While in these experiments, with one administration of cocaine, we anticipate the neuronal network to revert back to its initial state.With additional, consistent cocaine administration over time, we anticipate the network being driven to a new stable state becomes permanent.While experimenter-administered cocaine has been shown to produce different behavioral and neurobiological changes in animals compared to self-administration paradigms, this work is a crucial step in understanding initial neuronal network changes in response to cocaine.
Finally, how do our results help bring quantification to the current theories of cocaine affecting the brain.As the medial prefrontal cortex inhibits risky behavior in a top-down manner, one may be concerned with how functionality of the medial prefrontal cortex is altered to be able to "apply the brakes," so to speak, on other brain regions affecting behavior.Given our mesoscale experiments and analysis, it is still not clear how the functionality of the medial prefrontal cortex is altered, as this neuronal network is embedded in a sea of other neurons.Shutting down, or inhibiting, the functionality of the medial prefrontal cortex seems to be the intuitive stance and supported by physiological measurements at the brain region scale (15,16).We do indeed observe a partial shutting down, if you will, with a decrease in the neural activity amongst this neuronal subnetwork within the prefrontal cortex.However, as we observe changes in both functionalexcitatory and functional-inhibitory functionality at the individual neuron scale to arrive at different neuronal firing patterns in the presence of stimulation from external, or hidden, neurons.In other words, how changes at this mesoscale affect neuronal functionality at the larger brain region scale and then between brain regions must be understood, not only for cocaine, but for other disease states as well.To address this need, we incorporated the modern Hopfield model with high quality in vivo single neuron calcium imaging in a freely behaving rat.The development of technical and computational tools for neuroscience, such as the one we have described here, contributes to our ability to identify therapeutics that are able to restore dysregulated neuronal network function and advance human health.

Fig. 1 .
Fig. 1.Schematic and data from the experiments.Top: Implanted microscope with image of neurons fluorescensing.Bottom: Representative relative calcium intensity traces time prior and post cocaine administration for several neurons; partially created from Biorender.com.

Fig. 2 .
Fig. 2. RC model of a single neuron.a) The battery represents the voltage from all other neurons as well as any external signal.The switch S represents charging and discharging cycles of the neuron.When a neuron consumes ATP to open the ion channels, S is connected to the battery.During the discharge phase S moves from the battery to the other side to form a closed self-circuit.In this phase regardless of the incoming signals from the other neurons, the neuron does not fire.b) The current and voltage cycle of a neuron after firing.
We separately perform the minimization in Eq. 18 for prior and post cocaine administration.The predictions of Eq. 16 are compared with their corresponding true values in Fig.S1.The figure indicates that the test error is small and our predictions are close to the true values in most of the times.Our estimations for T ij and V (ext.)i for before and after cocaine administration are the outputs of this minimization.In other words, the sum of inputs V (ext.)i

Fig. 3 .
Fig. 3.The functional connectivity matrix functional-inhibitory and functional-excitatory neurons prior a) and post b) cocaine administration.The T ij matrix is derived using a time series machine learning analysis of the calcium signaling in the neuronal network with negative weights representing functional-inhibitory connections and positive weights representing functional-excitatory connections.

Fig. 4 .
Fig. 4. Changes in T ij to highlight the differences in the functional connectivity prior and post cocaine administration.a) The histogram demonstrates that there both increases and decreases in the weight connections.b) The changes in T ij prior and post cocaine.

Fig. 5 .
Fig. 5. Perturbing the neuronal network.a) A time-varying external current is applied to all neurons in the neuronal network prior to cocaine administration.b) The time-varying current of all the neurons of the network induced by the same external signal.Tsame external signal induce different currents in the neurons given the functional connectivity weights.Moreover, when the external current is disconnected, all induced current begin to decay, each with their own decay rate.

Fig. 7 .
Fig. 7.The inferred functional-excitatory and functional-inhibitory neuronal networks prior and post cocaine administration.a) and b) Networks represent prior to cocaine administration.c) and d) networks represent post cocaine administration.a) and c) Networks represent the functional-excitatory connections, while b) and d) networks represent the functional-inhibitory connections.The colors define the modularity classes, or clusters of closely interconnected nodes, using a community-finding algorithm outlined in Ref. (59) and implemented in Gephi(60).Moreover, the node sizes show the betweenness centralities.The interaction strengths are represented by the width of the edges.

Table 1 .
The betweeness centrality scores for the largest top five changes in neurons for both functional-excitatory (first five rows) and functional-inhibitory (second five rows) contributions.

Table 2 .
The learned parameters for before and after cocaine administration.