Modeling Synaptic Effects of Anesthesia and its Cortical Cholinergic Reversal

General anesthetics work through a variety of molecular mechanisms while resulting in the common end point of sedation and loss of consciousness. Generally, the administration of common inhalation anesthetics induces decreases in synaptic excitation while promoting synaptic inhibition. Animal studies have shown that, during anesthesia, exogenously induced increases in acetylcholine-mediated effects in the brain can elicit wakeful-like behavior despite the continued presence of the anesthetic. Less investigated, however, is the question of whether the brain’s electrophysiological activity is also restored to pre-anesthetic levels and quality by such interventions. Here we apply a computational model of a network composed of excitatory and inhibitory neurons to simulate the network effects of changes in synaptic inhibition and excitation due to anesthesia and its reversal by muscarinic receptor-mediated cholinergic effects. We use a differential evolution algorithm to fit model parameters to match measures of spiking activity, neuronal connectivity, and network dynamics recorded in the visual cortex of rodents during anesthesia with desflurane in vivo. We find that facilitating muscarinic receptor effects of acetylcholine on top of anesthetic-induced synaptic changes predicts reversal of the neurons’ spiking activity, functional connectivity, as well as pairwise and population interactions. Thus, our model results predict a possible neuronal mechanism for the induced reversal of the effects of anesthesia on post synaptic potentials, consistent with experimental behavioral observations. Author Summary Here, we apply a computational model of a network composed of excitatory and inhibitory neurons to simulate the network effects of changes in synaptic inhibition and excitation due to anesthesia and we investigate the possibility of its reversal by muscarinic receptor-mediated cholinergic effects. Specifically, we use a differential evolution algorithm to fit model parameters to match dynamics recorded in the visual cortex of rodents during anesthesia with desflurane in vivo. We find that changes of the fitted synaptic parameters in response to the increasing desflurane concentration matched those established by neurophysiology. Further, our results demonstrate that the cellular effects induced by anesthesia can be mitigated by the changes in cellular excitability due to acetylcholine.

exogenously induced increases in acetylcholine-mediated effects in the brain can elicit 23 wakeful-like behavior despite the continued presence of the anesthetic. Less investigated, 24 however, is the question of whether the brain's electrophysiological activity is also 25 restored to pre-anesthetic levels and quality by such interventions. Here we apply a 26 computational model of a network composed of excitatory and inhibitory neurons to 27 simulate the network effects of changes in synaptic inhibition and excitation due to 28 anesthesia and its reversal by muscarinic receptor-mediated cholinergic effects. We use 29 a differential evolution algorithm to fit model parameters to match measures of spiking 30 activity, neuronal connectivity, and network dynamics recorded in the visual cortex of 31 rodents during anesthesia with desflurane in vivo. We find that facilitating muscarinic A4/B1-B4 denote optimal parameter sets fit to experimental recordings at varying 221 anesthetic concentrations (0%, 2%, 4%, 6% desflurane, respectively reversal results. Parameters from initial fit used to simulate anesthetic effects and 229 cholinergic reversal. A/B-series describe optimal values of initial fit determined by the 230 differential evolution algorithm for network connectivity parameters obtained when ACh 231 effects are assumed constant (i.e., gKs is constant; A-Series) and when ACh effects are 232 allowed to change with anesthetic concentration (B-Series). Px denotes scaled changes 233 in synaptic conductance's mediated by the x receptor (x = NMDA, GABA and AMPA) as 234 described in Table 1 We first characterized the changes in the mean neuronal spike frequency as well as the 257 shape of the neuronal spike frequency distributions as a function of anesthetic level in the 258 optimized model networks (Fig 3 and Fig 4A). We observed that the neurons generally 259 fired less, in both experimental data and the simulations, as a function of anesthetic 260 concentration. Also, the spread of neuronal firing frequencies decreased significantly with 261 increased anesthetic level, with the loss of the right skew observed in the wake cases 262 (0%, A1 and B1). Spike frequency decreased as a function of desflurane levels for both 263 parameter series, (A and B series, without and with ACh changes, respectively), with a 264 similar frequency drop, irrespective of the implemented ACh changes that affect neuronal 265 excitability in the B series. In predicted ACh-induced reversal, the rightward skew in 266 frequency distributions was recovered, and the B series showed stronger recovery in 267 mean spike frequency as compared to the A series. This is because, as mentioned above, 268 accounting for cholinergic changes on neuronal excitability under desflurane anesthesia 269 predicts that synaptic changes are less severe. Namely, in the B series, GABAA synaptic 270 strength was not as high, and NMDA synaptic was not as low compared to the A series. Changes in experimentally recorded firing rate distributions under increasing desflurane 275 concentration (0, 2, 4, and 6%) show increased right skewness for the awake state in 276 comparison to anesthetic states. The bins were normalized by the total number of spikes 277 relative to the awake case (0% In both experimental and simulated results, the common feature was an increase in 287 network synchronization as a function of increased desflurane levels. Mean phase 288 coherence (MPC) measures the consistency of the relative phase that neurons fire with 289 respect to each other thus taking into account non-zero time lag synchrony. For both the 290 was lowest at the 0% anesthesia case and increased with increasing anesthetic 292 concentration ( Fig 4B). For higher levels of anesthesia, while the simulations and the 293 experimental data showed increasing MPC, simulated networks exhibited a larger 294 increase. This could be due to the fact that our model only represents local network 295 interactions, without incorporating the existence of external inputs that could additionally 296 desynchronize the network activity. In the visual cortex, there are non-local network 297 inputs possibly preventing a high level of synchronization in the locally recorded network 298 activity. 299

300
The anesthetic reversal with increased levels of ACh (i.e. decreased gKs) led to decreases 301 in MPC indicating desynchronization in network activity. This is not suprising, as 302 decreased M-current leads to larger differentiation in neuronal firing rate as a function of 303 external input and changes in the phase response curves (PRCs), which also promote 304 desynchronisation for lower gKS [21]. 305 306 Anesthetic and reversal effects on network information metrics 307 We computed the information theoretic measures network integration (I(X)) and 308 complexity (C(X)) for both experimental data and simulated network activity. Integration 309 I(X) is a generalization of mutual information that measures the amount of total entropy 310 of a system that is accounted for by the interactions among its elements. I(X) is zero when 311 system elements are statistically independent [22]. Complexity C(X), on the other hand, 312 measures the total entropy loss due to interaction of system elements, or, equivalently, of the entire system. C(X) is low for systems with independent elements or with highly 315 synchronous elements.  We estimated network excitatory and inhibitory synaptic strengths, as well as network 359 excitatory and inhibitory connection probabilities, in the optimized networks and 360 compared them directly to these same measures computed from the experimental data. The results discussed above report trends observed for measures of average network 408 activity, such as frequency, mean phase coherence, integration and complexity, as well 409 as network connectivity strength and probability. And while ACh-induced reversal 410 reinstated these network-level measures, the measures do not account for recovery of 411 functional connectivity in the network which would contribute to information processing. 412 In this section, we investigate how ACh reversal affects the relative frequency profile of 413 individual neurons with respect to other neurons in the network and also look at effects of whether the internal dynamic structure of network activity is reinstated during the ACh 416 reversal. 417 experimental data and in the optimized networks at each level of anesthetic concentration 420 to its firing rate in the waking state (Fig 6 and S1 Fig). In the figure panels, the x-axis 421 represents firing frequency of individual cells for different anesthetic levels and the y-axis 422 represents the firing frequency for the same cells in the non-anesthetic (0% or A1/B1) 423 conditions. For the experimental data mutliple units can be potentially detected on a single 424 electrode. This led to potential ambiguity in neurons assigned across anesthetic levels. 425 To address the, neuron identity was based on firing rate in the 0% case. Namely, for units 426 We observed that, generally, in both experiments and simulation results the relative 434 frequency of the neurons was preserved, i.e. neurons that fired at higher frequencies as 435 compared to other cells in non-anesthetic conditions retained higher firing frequencies at 436 the different anesthetic levels, albeit absolute frequencies decreased. Conversely, anesthetic state continued firing at lower relative frequencies in the anesthetic conditions. 439 Qualitatively similar results were observed for A-series networks (Fig 6) and B-series 440 networks (S1 Fig). In our study we found that potentiation of inhibitory GABAergic and inhibition of excitatory 529 glutamatergic NMDA synaptic receptors do indeed lead to graded decreases in 530 population activity and increases in synchronization, as quantified by firing rate and mean 531 phase coherence, as well as measured decreases in integration and complexity. 532 Additionally, we were able to recover changes in functional network connectivity which 533 matched changes seen in literature [25,26]. The simulation results were robust; although 534 some of the measures (frequency, MPC, I(X) and C(X)) were used for optimization of 535 model parameters via the differential evolution algorithm, the results held for a wide range 536 of non-fitted measures within physiologically reasonable limits. Moreover, the parameter 537 fits obtained for increasing levels of anesthetic matched in their relative magnitudes to the 538 reported anesthetic induced changes in synaptic efficacy. 539

540
Understanding the mechanism of anesthesia through computational modeling 541 The cellular mechanism of anesthetic action with respect to loss of awareness has been 542 a subject of intense investigation. Computational models are actively used to make 543 progress in this area of research. Because differing classes of anesthetics elicit different 544 effects on synaptic receptor subtypes, many modeling approaches aim to determine how 545 nuanced changes in receptor binding and synaptic activity lead to changes in neural or glutamatergic synaptic changes are attributed to a single parameter that maps to different 548 concentrations of general anesthesia [27]. Other modelling approaches seek to 549 understand the mechanism of specific anesthetic agents; for example, the effects of 550 propofol have been studied through the modeling of both GABAA and GABAB 551 amplitude/duration and the effects on cortical synchrony and EEG rhythms [18,28]. 552 Enflurane and isoflurane are other commonly modeled anesthetics where the roles of 553 both glutamatergic receptor binding and GABAergic effects are taken into consideration 554 [28][29][30]. Anesthetic action effected through post synaptic potential (PSP) changes, from 555 a modelling perspective, is a relativity robust explanation supported by its effectiveness 556 across modelling paradigms. These include "mean field" models as well as networks of 557 "integrate and fire", "Izhikevich "and "HH" neurons, which all show reduced activity and 558 changes to population synchrony when modeling anesthetic effects on synaptic receptors 559 [30 -32]. 560

561
Our study is distinguished from former computational models of anesthetic effects by the 562 independent consideration of the effects on NMDAR and GABAR through PSP changes, 563 as well as of cholinergic influence through changes in the muscarinic M-current. We also 564 used a more biologically realistic log-normal distribution for synaptic weights [33]. 565 Because we had access to experimental spike data, we were able to directly fit our model 566 to empirical data at graded levels of anesthesia and then test our hypothesis regarding 567 cholinergic anesthesia reversal. GABAergic pathways (Fig 1). By considering solely the effect of changes on IPSPs via 631 GABAR and EPSPs through NMDAR we show that not only can changes in population 632 activity (firing rate, synchronization and entropy), be accomplished without changes in 633 cholinergic influence but that increasing cholinergic influence alone can reverse these 634 effects. This demonstrates that cortical cholinergic presence has the potential to mitigate 635 the general effects of inhalational anesthesia. In many cases, however, such as for the 636 effects of desflurane, inhalational anesthesia can affect muscarinic and nicotinic ACh 637 receptor binding and for this reason we decided to model the cooperative effects from case of cholinergic reversal, however, this confounded the role of ACh, as the changes 640 in ACh due to anesthesia could be argued to be trivially reversed in the reversal states. 641

642
In this study, we used measures of synaptic functional connectivity, computed from 643 average pairwise correlations of neuron spiking, to quantify changes in overall network 644 behavior in both anesthesia and reversal conditions. We showed that the cosine similarity 645 in the functional connectivity matrix increased for the full reversal state when compared 646 to the high anesthetic state. This means that specific neuron to neuron functional 647 connectivity was highly correlated between the awake and reversal states but not the 648 anesthesia states. This suggests that the functional topology of a network can be 649 reversed through a different receptor pathway than is used to achieve the state of [15], and a similar reversal from isoflurane was observed in response to microinjection of 660 histamine into the basal forebrain [49]. Unlike general anesthesia, however, the 661 gabazine, on the effects of propofol as well as ketamine [50]. The application of gabazine 664 led to wake-like responses when rats were sedated with propofol, which acts through 665 potentiation of GABAA receptors, but gabazine was ineffective when used during 666 administration of ketamine, which has been known to act through modulation of NMDA In summary, we demonstrated that experimentally observed changes in neural activity 714 and functional connectivity caused by desflurane could be computationally reproduced by 715 modulating synaptic efficacy according to the known synaptic effects of the anesthetic. 716 Additionally, we showed that by modulating the M-current alone, the effect of anesthesia 717 on neural activity and functional connectivity in the network could be at least partially 718 reversed. In the future, more comprehensive models that take into account cortical 719 architecture, thalamocortical interactions and a broader array of cellular mechanisms will 720 help to fully understand the complex roles of synaptic modulation in producing the 721 observed neuronal network and behavioral effects of anesthesia. Neurons were connected randomly with 10% probability. Synaptic strengths followed a 783 log normal distribution, as suggested to occur in cortical networks (Fig 8A) [33] .The information that measures the amount of total entropy of a system that is accounted for 832 by the interactions among its elements. I(X) is zero when system elements are statistically 833 independent [22]. C(X) measures the total entropy loss due to interaction of system 834 elements, or, equivalently, the difference between the sum of the entropies of the 835 individual elements and the entropy of the entire system. C(X) is low for systems with 836 independent elements or with highly synchronous elements. 837 To compute these measures, the total spiking activity from an experimental recording or 838 a network simulation was partitioned into patterns by binning spike trains into 1 ms time 839 bins and constructing vectors for each time bin containing a 1 at the neuron index if the 840 neuron spiked within that time bin and a 0 if there was no spike (columns in Fig 9). The 841 set X of unique vectors, representing patterns of spiking activity within a bin, that occurred centered on each spike of the designated "reference" cell of the pair and discretized into 893 1.3 ms bins. Cross-correlations of discretized segments between the "reference" and 894 "comparison" cell for every "reference" cell spike were summed to form cross 895 correlograms (Fig 10). between example pairs of "reference" and "comparison" cells, centered at spike times of 904 the "reference" cell, from the experimental recordings (left column) and simulated 905 networks (right column). Significance bands were computed from a jittered data set of 906 "comparison" cell spike times (gray line = mean of jittered data set, red line = excitatory 907 significance, blue line = inhibitory significance, see text). A-B) Example cross 908 correlograms showing significant excitatory connections between cell pairs. C, D) 909 Example cross correlograms showing significant inhibitory connections between cell 910 pairs. 911 912 913 repeated by 100 times to for the jittered data set. The global confidence band for 917 excitatory (inhibitory) connectivity was computed by taking the 97% confidence interval 918 associated with the global peak (trough) of the jittered data set. A significant connection 919 was determined when the peak (trough) of the original cross correlogram was greater 920 (less) than 2 times the 97% confidence interval when measured from the mean (blue/red) 921 The 20 lowest cost parameter sets were kept and each parameter value was randomly 953 varied uniformly by 10% of its value to avoid duplicate values. The final 10 parameter sets 954 were then constructed using the differential evolution algorithm. 955 956 Similar to typical differential evolution procedures [66,67] we set a cross over probability 957 CR = 0.8 and a differential weight DW = 2. From the subpopulation of 20 parameter sets, 958 ( ( = 1, … ,10). For each set ( , 3 different sets ( , ( and were chosen that were 960 different from ( and each other. Then, for each element = 1, … ,4 in the set, a random 961 numberfrom the uniform distribution [0,1] was chosen. Ifwas less than CR, a new 962 parameter value -( was generated as - This was done for each element in the new agent and was repeated until 10 new agents 964 were created. After this was done the 10 new agents were simulated and then the 30 total 965 parameters were evaluated for their cost. The 10 with the highest cost (worse fit) were 966 then rejected and the process was repeated. uniformly by +/-5% to generate variability in the event of parameter convergence. In the 979 optimizations for the 2% and 4% anesthetic cases, A2/B2 and A3/B3, respectively, the 980 initial population for A2/B2 was A1/B1, and the initial population for A3/B3 was A4/B4. Evolutionary algorithm procedure, differential evolution, was used to optimize model 985 parameters. For each generation, 10 agents (parameter sets) with the highest cost 986 function from the population of 30, were chosen for replacement. Algorithm was repeated 987 until stopping criteria of 100 generations without change in lowest cost function value 988 across the population was met. B,C) Lowest cost function values across the parameter 989 set populations at each generation for the A-series (B) and B-series (C) parameter 990 optimizations. Population A1/B1, A2/B2, A3/B3 and A4/B4 were optimized to 991 experimental data from the 0%, 2% ,4% and 6% anesthetic cases, respectively. The 992 optimizations for A1 and B1 were identical. In the A-series (A2-A4), PNDMA, PGABA, were 993 optimized and in the B-Series (B2-B4), PNDMA, PGABA, gKs were varied 994 995 . 996 network realizations, keeping the network structure fixed for all anesthetic levels. The 999 average and error (SEM) for the optimized parameters across these 10 networks is shown 1000 in Table 1. Table 2 lists the parameter values with the lowest cost function for one of these 1001 optimization runs that we used in our model analysis.