In-silico evidence that increased cortico-thalamic connectivity and superficial pyramidal disinhibition underlie broadband task-related spectral changes induced by ketamine

Cortical recordings of task-induced oscillations following subanaesthetic ketamine administration demonstrate alterations in amplitude, including increases at high-frequencies (gamma) and reductions at low frequencies (theta, alpha). To investigate the population-level interactions underlying these changes, we implemented a thalamo-cortical model capable of recapitulating broadband spectral responses. Compared with placebo, ketamine was found to increase the decay rate of NMDA receptor currents, decrease the decay rate of GABA-B receptor currents, but had no effect on AMPA or GABA-A currents. Furthermore, ketamine decreased the inhibitory self-modulation of superficial pyramidal populations, increased the inhibitory self-modulation of inhibitory interneurons and increased the strength of the cortico-thalamic projection from layer 6 into thalamus. In-silico rectification of each of these parameters to their placebo state revealed the accompanying spectral changes. Similar broadband effects were demonstrated for the NMDA constant and intrinsic connectivity changes, whereby rectification resulted in normalising the amplitude of alpha (increased) and gamma (decreased). While supporting theories of superficial pyramidal disinhibition following ketamine administration, our results suggest a role for altered cortico-thalamic connectivity relevant to understanding ketamine induced cortical responses and potentially identifies a system-level mechanism contributing to its antidepressant effects.


Introduction
The spectral composition of magnetoencephalographic (MEG) signals derived from visual cortex reflect a complex system of neuronal interactions, revealed as frequency specific oscillations 1,2 and power law-like phenomena. Individual differences in common metrics, such as the peak frequency and amplitude within specific frequency windows, can predict behaviours 3 , disease states [4][5][6][7][8][9] , and are sensitive to experimental manipulations, including pharmacological 10-15 .
While the laminar and network generators of frequency-specific oscillations remain equivocal, prevailing theories propose distinct dominant laminar generators for oscillations of different frequencies. One particular proposal is that high frequencies tend to be formed by current densities in supragranular layers, while lower frequencies appear to predominate in deeper regions of cortex, around layer V/VI 20,25 or in thalamo-cortical loops 26,27 . Under this model, the noted effects of ketamine on both low and high frequency bands suggest altered neuronal functioning across much of the cortical laminae, and possibly in thalamo-cortical connectivity.
Computational modelling is an approach based on the use of basic mathematical descriptions of neuronal units, such as neural mass models, to explore neurophysiological phenomena. Methods such as Dynamic Causal Modelling (DCM) [28][29][30][31] provide frameworks for generating typical and atypical empirical features in imaging and electrophysiological time series. In DCM the subsequent fitting (inverting) of parameterised mean-field neuronal models to empirical data features, acquired under different pharmacological or task states, permits an in-silico assay of usually unobservable neuronal states such as synaptic connectivity between cell populations, or the decay times of specific receptor types. This approach has received construct and face validations 32,33 .
We implemented a conductance-based thalamo-cortical circuit with realistic intrinsic dynamics, to assay both changes in cortical and thalamo-cortical population connectivity and receptor dynamics underlying observed spectral changes induced by a subanaesthetic, anti-depressant relevant dose of ketamine, as recorded with magnetoencephalography (MEG). We hypothesised that alterations in NMDA receptor activity would be identified under this manipulation and further sought to identify additional effects that might underlie Ketamine's complex dose-dependent role as a psychotomimetic and anti-depressant. cortex at the location of peak t-statistic. Spectra were computed using a fast Fourier transform of the virtual sensor.

Ketamine and placebo infusion protocol
All subjects underwent both a ketamine and a placebo infusion, in counterbalanced, pseudo-random order, with at least 2 weeks between sessions. The ketamine infusion consisted of 0.5mg/kg (of body mass) racemic ketamine hydrochloride in 50 ml saline. For placebo, 50 ml saline only was used.
Infusions were administered by intravenous cannula in the dorsum of the left hand, with an initial bolus (~ 1 minute) of 0.25 mg/kg followed by maintenance infusion of .25 mg/kg delivered over 40 minutes.
Maintenance infusion was controlled by an Asena-PK infusion pump (Alaris Medical, UK).
Populations were connected in a laminar structure in accord with their anatomical position and connections though the model does not explicitly encode spatial distribution. This meant we had pyramidal and interneuron populations in layer 2/3 and layer 5, a stellate population in layer 4, and a thalamic-projection pyramidal population in layer 6. The time evolution of each population was governed by coupled differential equations describing the rate of change of membrane potential and rate of change of conductances supported by particular receptors; Where V is the membrane potential, ) is the conductance from a particular receptor type and ) ( − ) ) represents the contribution to the membrane potential of channel k, C is the membrane capacitance, ) reversal potential of channel k and u is input current.
In the conductance equation, ) is the decay rate of currents mediated by channel k and ) represents presynaptic input (firing rate x connectivity -see below). Conductance terms are included for 6 channels, AMPA, NMDA, GABAA, GABAB, m-currents and h-currents, from voltage-gated and leak potassium channels respectively. Only layer 6 thalamic-projection pyramids and thalamic relay populations possess m and h channels (see figure 1). This formulation is an extension of the NMDApossessing model of Moran 39,41 (CMM_NMDA in the SPM DCM toolbox) and includes the same voltage switch to represent the NMDA magnesium block; Reversal potentials and decay rates (table 1) were obtained from existing models 39,41 and from the literature (GABAB 42 , M-channel 43,44 , H-channel 45 ). M-channels are responsible for m-currents; noninactivating voltage-gated potassium currents, which may be important to the intrinsic dynamics (oscillations) of the cell membrane 44 . H-channels produce hyperpolarisation activated currents, which not only have a role in the regulation of the resting membrane potential 46 but may contribute to macroscopic network oscillations 45 . The G-protein coupled GABAB receptor has also been linked to macroscale oscillations, including in the gamma range 47 , possibly through its ability to re-balance the excitation-inhibition balance through interactions with NMDA receptor mediated currents 48 . For each population in the model, the input to population i, is the expected firing rate from the source population j; the coupling of this presynaptic input to a postsynaptic response is represented by the parameter γ(i,j) which subsumes various biophysical processes such as receptor binding and transmitter reuptake.
Here, the sigmoid function sigma represents the cumulative distribution function of the presynaptic depolarisation around a threshold potential of VR = -40 mV, which determines the proportion of cells firing 39,49 .
The model was integrated numerically over a 3 s period (ensuring steady-state) at dt = 1/1200 s using an Euler method with delays as per David et al 50 . Input was a constant direct-current to thalamic relay cells. The cortex-thalamus-cortex loop took 80 ms, made up of 60 ms cortex to thalamus and 20 ms thalamus to cortex, in accordance with previous models 51, 52 . The spectral response was calculated by Fourier transform of a weighted sum of the population membrane potentials (superficial pyramidal * 0.8 + spiny stellate * 0.2 + layer 5 pyramidal * 0.2 + layer 6 pyramidal * 0.2). The first 300 ms of simulation was disregarded.
The model was fitted to the MEG virtual sensor frequency spectrum between 4 and 80 Hz using the DCM Gauss-Newton optimisation routine (a Variational Laplace / Bayesian inversion routine for nonlinear models, see spm_nlsi_GN.m), which minimises a free energy term [53][54][55] in order to maximise the model evidence and produce a-posteriori estimates for the model parameters, which included intrinsic connectivity parameters ( in eq.3) and receptor time constants for AMPA, GABAA, NMDA and GABAB ( in eq.1).

Study objectives
In this study, we had 2 key objectives. First, to assess how well the model recapitulated the broadband spectrum of the MEG virtual sensors, and whether the model fits were adequately accurate to recapitulate the drug effects on the spectrum.
Second, given adequately fitted models, we were specifically interested in examining the effects of ketamine, as compared to placebo, on 2 sets of model parameters -as well as how these parameter effects change the output spectrum of the model. The parameter sets included (1.) the time-constants of the modelled receptors (AMPA, NMDA, GABAA and GABAB; not m-or h-channels as there were not parameterised but reflect intrinsic population dynamics) and (2.) the synaptic connectivity between populations, as reflected in figure 1(A). Between drug differences were determined by paired-t test with 5000 permutations and omnibus correction for multiple comparisons 56 . For the time-constants, this correction was for the 4 TC parameters. For the synaptic connectivity parameters, this correction was for 22 parameters (those depicted in figure 1(A), minus the 2 connections originating from reticular populations, which were fixed because they had no effect on the model output).

Participants
Of 20 recruited subjects (mean age 25.7, SD 6.2), 1 subject withdrew prior to ketamine infusion, and an error in MEG acquisition was made for another. A further 2 datasets from the visual paradigm reported here, were rejected based on visual inspection of gross recording artifacts (reported previously 57 ), resulting in 16 usable, complete (i.e. including ketamine + placebo) datasets for analysis.   Three connectivity parameters were different between ketamine and placebo conditions; inhibitory SPàSP self-modulation decreased with ketamine (t = -3.01, p = 0.020), inhibitory SIàSI selfmodulation increased with ketamine (t = 2.75, p = 0.024) and the cortico-thalamic projection from layer 6 pyramids into thalamic relay populations was increased with ketamine (t = 3.13, p = 0.012). Performing the same quantification for the rectified GABAB TC revealed the same pattern of changes, but with small effect (mean alpha amplitude +0.25% at 10 Hz, mean gamma -0.02% peaking at 48 Hz, see figure 7B).

Discussion
Model sensitivity In this study, we have demonstrated that a simplified thalamo-cortical model, implemented with DCM, is able to recapitulate the spectral changes induced by subanaesthetic ketamine administration during a visual grating task. Comparing the models fitted to ketamine and placebo, we identify 5 key parameters that were changed by ketamine; namely, the receptor decay constants of NMDA (shortened) and GABAB (lengthened) channels, the inhibitory self-modulation of superficial pyramidal (decreased) and inhibitory interneurons (increased) and the cortico-thalamic projection from layer 6 into thalamic relay populations (increased).
We assessed the goodness-of-fit of the model spectra using correlations over all datasets of the amplitude at each frequency step, comparing the model and 'real', MEG virtual sensor spectra. As evidenced in figure 3, the model accurately summarized the features of the real spectrum, with significant correlations through much of the frequency spectrum, declining only at higher frequencies (> 70 Hz) and only non-significantly at 80 Hz.
Before making any between-drug parameter inference, it was important to assess whether the model had summarized the data well enough to be sensitive to the drug induced changes in the spectrum. As shown in figure 4, the spectral difference in amplitude observed in the empirical MEG data (figure 4A), which included significant reductions in alpha amplitude and increases in gamma amplitude, were also evident at the same frequencies in the model output ( figure 4B).

Ketamine parameter effects -NMDA channel decay time constant
While its pharmacological profile is known to be complex 58  processing and alters gamma oscillations 70,72,75 . It has been proposed as a candidate mechanism by which ketamine disrupts gamma 71 , a theory which our result could support. -

Increased inhibition within SI populations (SIàSI)
Ketamine was found to increase the inhibitory self-modulation of superficial interneuron populations, compared to placebo. Initially, this result seems to contrast the theory of disinhibition outlined above, because the activity of inhibitory interneurons should be decreased 70 . This contrast may reflect the simplistic parameterisation of cortical models like ours, where every population has its own inhibitory self-modulation parameter as well as their being explicit inhibitory interneuron populations. Under this framework, increased inhibitory self-modulation of inhibitory interneurons (SIàSI) is an unrelated process to the inhibitory self-modulation of superficial pyramidal populations (SPàSP).
-Increased cortico-thalamic output connectivity (TPàRL) Ketamine increased the strength of the excitatory cortico-thalamic projection from the layer 6 pyramidal population into thalamic relay population, compared to placebo. Ketamine mediated increases in cortico-thalamic connectivity have been reported in resting state functional magnetic resonance imaging (fMRI) studies 76 , although this study focussed on somatosensory and temporal cortical projections 76 , rather than visual.
In rat, administration of ketamine was found to enhance the visual cortex layer 6 projection to thalamus 77 , mimicking exactly the alteration found in our model. They also demonstrated that this connection predicted the amplitude of thalamic gamma oscillations 77 , however qualitative assessment of the output of our fitted models did not reveal oscillations in the gamma range in the membrane potential of thalamic populations.

Parameter rectification and spectral effects
For each of the previously described 5 parameter changes, we quantified the effect that restoring that parameter (to its value from the placebo model) had on the output spectrum. This in-silico 'reverse lesioning' allowed us to understand the effect that each of the ketamine-altered parameters had on the output of the model. Or, which parameter changes might underlie which spectral changes induced by ketamine.
Quantifying the parameter changes as percent-change from the baseline fitted ketamine model c.f. a sensitivity analysis, all five parameters demonstrated similar effects of increasing the amplitude of low frequency oscillations (which had been reduced by ketamine) and reducing the amplitude of high frequency, gamma oscillations (which had been increased by ketamine). Effects of the NMDA and GABAB TCs were to increase the 10 Hz alpha by +1.6% and +0.25%, respectively. The NMDA TC reduced gamma with peak effect at 50 Hz, by -0.2%, while the GABAB TC reduced gamma with a peak effect at 48 Hz, by 0.02%. The 3 connectivity parameters all also increased alpha at 10 Hz, but with varying peak frequency effects (SPàSP 45 Hz -0.02%, TPàRL 51 Hz -0.07% and SIàSI 55 Hz, -0.12%).

Study limitations and future work
Given the sensitivity of model parameters to manipulation by ketamine, it would be interesting to know whether any of the observed effect are dose-dependent, particularly since ketamine has seen a resurgence in popularity due to its antidepressant effects. In the present study, we used a fixed 0.5mg/kg dose, however there have been no neuroimaging-based studies of ketamine at varying levels, which would permit discovery of dose-dependent effects.
This study explored the effects of ketamine in healthy, young participants. However, the same modelling framework could be applied to clinical trial data, where ketamine is currently under investigation as a treatment option for a variety of psychiatric (mood) disorders. In this context, the sensitivity of model parameters could serve as sensitive and specific biomarkers -both of illness and of treatment response.
In this study, we have demonstrated that a thalamo-cortical circuit, fitted to empirical pharmaco-MEG data with DCM, is accurate enough to recapitulate spectral changes induced by ketamine; and sensitive to local receptor and synaptic connectivity dynamics. We have shown that both superficial layer disinhibition and cortico-thalamic afferents play a role in shaping the oscillatory responses of ketamine, which should be considered in future studies linking gamma and disinhibition.