A multivariate physiological model of vagus nerve signalling during metabolic challenges in anaesthetised rats for diabetes treatment

Objective. This study aims to develop a comprehensive decoding framework to create a multivariate physiological model of vagus nerve transmission that reveals the complex interactions between the nervous and metabolic systems. Approach. Vagus nerve activity was recorded in female Sprague-Dawley rats using gold hook microwires implanted around the left cervical vagus nerve. The rats were divided into three experimental cohorts (intact nerve, ligation nerve for recording afferent activation, and ligation for recording efferent activation) and metabolic challenges were administered to change glucose levels while recording the nerve activity. The decoding methodology involved various techniques, including continuous wavelet transformation, extraction of breathing rate (BR), and correlation of neural metrics with physiological signals. Main results. Decrease in glucose level was consistently negatively correlated with an increase in the firing activity of the intact vagus nerve that was found to be conveyed by both afferent and efferent pathways, with the afferent response being more similar to the one on the intact nerve. A larger variability was observed in the sensory and motor responses to hyperglycaemia. A novel strategy to extract the BR over time based on inter-burst-interval is also presented. The vagus afferent was found to encode breathing information through amplitude and firing rate modulation. Modulations of the signal amplitude were also observed due to changes in heart rate in the intact and efferent recordings, highlighting the parasympathetic control of the heart. Significance. The analytical framework presented in this study provides an integrative understanding that considers the relationship between metabolic, cardiac, and breathing signals and contributes to the development of a multivariable physiological model for the transmission of vagus nerve signals. This work progresses toward the development of closed-loop neuro-metabolic therapeutic systems for diabetes.


Introduction
Glucose control is a challenging task in diabetes because a complex network of mechanisms including hormones, metabolites, and neurotransmitters interact with each other to ensure a robust and tight regulation of blood glucose fluctuations.The traditional therapeutic strategy in type 1 diabetes involves purely metabolic control.The current state of the art technology, called the artificial pancreas, replicates the endocrine function of a healthy pancreas and calculates the optimal dose to be delivered based on the real-time glucose readings from a sensor.Despite the improved glycaemic control recently reported in clinical trials (Brown et al 2019), the system deviates from the physiological control insofar as it has no anticipatory and fast control components.This results in delays in insulin action and risks of hypo and hyperglycaemia, which correspond to low and high-glucose levels respectively (Ramli et al 2019).
Recent advances in multidisciplinary research have enabled the emergence of bioelectronic medicine, which uses electrical interfacing with peripheral nerves and promises to be a powerful therapeutic approach to overcome current challenges in the treatment of diabetes (Güemes Gonzalez et al 2020).In this context, the vagus nerve provides an ideal target to obtain additional metabolic information as it plays a major role in the regulation of multiple physiological processes by providing an integrative response coordinating and balancing the activity of many internal organs (Câmara andGriessenauer 2015, Karemaker 2022).This integrative response has been demonstrated in the last decades through the important role of the vagus nerve in regulating the immune system and coordinating the response to infections and inflammation (Pavlov and Tracey 2005, Joseph et al 2019, Pavlov et al 2020).
The specific contribution of the vagus nerve to control glucose levels is complex and multi-faceted, as it has a key role in regulating various aspects of glucose metabolism and hormone secretion.The vagus nerve transmits two key states of blood glucose levels: hyperglycaemia and hypoglycaemia.The vagus nerve transmits hyperglycaemia information through afferent glucose-inhibited neurons from a complex system of glucose-sensing cells located in different parts of the body to different brain regions including the amygdala, hindbrain and hypothalamus (Routh 2010, Grill and Hayes 2012, Güemes and Georgiou 2018, Belsham and Dalvi 2021).This, in turn, triggers efferent responses to control a variety of internal processes, such as glucose production and storage, fat breakdown and uptake in the liver, insulin release during early phases and optimising post-meal insulin secretion (Pocai et al 2005b, Thorens 2011), and regulation of glucose utilisation by peripheral tissues like muscle, brown and white fat (Lin et al 2021).
The vagus nerve plays an especially important role during hypoglycaemia, which is a risky situation that can lead to symptoms such as dizziness, confusion, and unconsciousness.The motor vagal neurons drive the counterregulatory response to hypoglycaemia by increasing hepatic glucose production-initially via glycolysis and later by gluconeogenesis (Pocai et al 2005b)-, increasing glucagon secretion from the pancreas (Lamy et al 2014), and impacting the gastrointestinal function (Hjelland et al 2005, Zhou et al 2008).Overall, all these processes work together to ensure that glucose levels are regulated within a narrow range, which is important for maintaining energy levels and supporting the functioning of other physiological systems in the body, including the cardiac, respiratory and digestive systems.
The intricate interplay of various physiological systems controlled by the vagus nerve makes it a formidable challenge to understand how it co-ordinately controls energy metabolism and homoeostasis.Therefore, the development of neuro-metabolic therapeutic systems that leverage the neural information conveyed by the vagus nerve requires an integrative understanding that considers the relationship between metabolic, cardiac and breathing signals.In this work, we examine if information about changes in each of these systems, and in particular glucoserelated clues, can be isolated from neural recordings from the surface of the cervical vagus nerve.We hereby present a framework to decode the motor and sensory vagal signalling of glucose, heart rate (HR) and breathing using descriptive data mining techniques, and contribute to the development of a multivariable physiological model for the transmission of vagus nerve signals.

Animal surgery
A total of 23 adult female Sprague-Dawley rats (Charles River Laboratories) weighing 200-300 g were used for acute electrophysiology experiments.In 11 animals the vagus nerve was kept intact, and vagotomy was conducted in the remaining animals for recording afferent activation from the distal end of the ligated nerve (n = 6) and for recording efferent activation from the proximal terminal (n = 6).Animals were received at least 1 week before surgery for acclimation and group-housed in individually ventilated cages with ad libitum access to food and water and a 12 h light cycle (7 am-7 pm).Food, but not water, was withdrawn 1.5 h prior to the start of the experiments.Animals were carefully handled using a tunnel to minimise the experience of stress.All animals were anaesthetised with an intraperitoneal (IP) injection of urethane (1.5 g kg −1 ) and were included in the analysis (no animal was excluded from the study).Urethane was selected as the anaesthetic agent as it does not depress neural responses when compared to volatile anaesthetic agents such as isoflurane (Silverman et al 2018).The temperature of the animal was monitored constantly and kept between 36 • C and 38 • C during the surgery and implantation of devices via a self-regulating heating pad.All work was carried out in accordance with the relevant ethical guidelines.
Exposure of the left vagus nerve was done following the surgical procedure described elsewhere (Güzel et al 2014).Animals were placed in a supine position and a 15 mm longitudinal cervical midline incision was made.The vagus nerve was exposed free from the carotid bundle by dissection of the perivascular structures (fascia, adventitia, and sympathetic plexus) using sharply curved forceps.Great care was taken in this step to avoid excessive manipulation or lesion of the nerve and to keep consistency in the degree of vagus nerve detachment across all animals.A region of 3 mm was exposed at the tip of polytetrafluoroethylene (PTFE)-coated gold microwires (75 µm diameter) and implanted around the epineurium of the nerve in the form of a 360 • hook.In animals undergoing ligation, the vagus nerve was held in place using a suture, the gold microwires were implanted and the nerve was cut distally or proximally for efferent and afferent recordings, respectively.A second electrode consisting of gold microwire was placed at one of several locations under the skin and on top of muscles close to the vagus nerve recording site, to capture and identify any non-neural signals, such as electromyography (EMG) artefacts.Once the surgery was completed, including nerve transection in the respective cohorts, a 15 min stabilisation period was allowed before starting the recording, which was enough for HR and body temperature to stabilise after surgery.Impedance values were always measured during the stabilisation period at 1 kHz using the Intan capabilities (in some cases the measurement was repeated at the end of the experiment with minimal variations).The recording period consisted of 30 min of baseline recording followed by control and metabolic challenges of 1 ml saline IP injection at minute 30 and 1 ml insulin (100 U ml −1 ) IP injection 20-50 min after the saline injection.In four animals, an IP injection of glucose (2 g kg −1 ) was delivered 40-50 min after the insulin injection (figure 1(a)).The total duration was approximately 120 min.Glucose levels were monitored every 5 min from puncturing superficial veins in the tail using a handheld glucometer (Abbot FreeStyle).This sampling rate was chosen based on a balance between capturing meaningful glucose dynamics and minimising the data collection burden over extended periods.

Acute electrophysiology
The recording devices were connected to a 32channel recording headstage (Intan Technologies, Los Angeles, CA, USA) via a custom-built omnetics/ZIF connector PCB and flexible flat cable.Signals were acquired at 30 kHz through the RHS 2000 stimulation/recording controller (Intan Technologies, Los Angeles, CA, USA).All the recording devices were grounded via a common gold film of 1 cm 2 coated by poly(3,4-ethylenedioxythiophene) polystyrene sulfonate (PEDOT:PSS) placed subcutaneously close to the recording site to allow for a single-ground point.This not only simplified the setup and reduced the number of connections needed, but also helped minimise the introduction of ground loops, which can cause electrical noise interference in the system.

Computational analysis
We have conducted a comprehensive multivariate descriptive analysis using advanced filtering and artefact reduction and data mining techniques of clustering and multivariable correlation analysis.

Cardiac component: wavelet decomposition for ECG detection
The cardiac signal was predominantly found in the recordings and therefore good filtering of the ECG signal was needed to avoid detecting cardiac beats as neural spikes (figure 1(c)).For this purpose, raw recordings were first low-pass filtered using a 4th order Butterworth (cut-off frequency at 100 Hz).The cardiac component of the resulting signal was isolated by applying continuous wavelet transformation using the Mexican hat wavelet family at a 5 ms scale, which resembles the duration of a heartbeat.The cardiac spikes were identified by applying a threshold that was calculated as: where σ ci is the standard deviation of the Gaussian noise of the cardiac signal ci, N is the number of samples and 0.6745 is the 75th percentile of the standard normal distribution (Donoho and Johnstone 1994).From the detected peaks, the R-R intervals (distance between R peaks in ms) were computed and consequently corrected by replacing any values with a z-score greater than 2 with the median of the array.
The HR evolution over time (in bpm) was eventually obtained from the corrected R-R signal (in seconds) using the following formula: HR (bpm) = 60/RR (s) (figure 1(e)).For a full overview of the processing of the cardiac signal see supp.figure 2.

EMG detection
High-amplitude spikes corresponding to EMG activity, such as fast muscle contractions or movement artefacts, were confirmed using a second microwire on top of the muscles surrounding the vagus nerve recording site.These spikes were identified by applying a threshold on the maximum amplitude for the detected spikes of five times the standard deviation (SD) from the mean and were subsequently discarded to not be considered as neural activity.

Inter-burst-interval and breathing rate
To evaluate the bursting activity in the raw neural data, we first computed the spectrogram from the power spectral density of the neural signal (frequency range [0-4000] Hz, NFFT = 2048).Then, for each time value, the contribution of frequencies greater than 500 Hz was summed, which clearly reflects multiple-unit spiking activity, within nonoverlapping windows of 1 s.The bursting activity was then clearly isolated from background activity and individual peaks were identified to determine the temporal distances between two following bursts (in seconds).A rolling window of 1 s duration (overlapping) was applied to extract the mean temporal distance, and the resulting signal comprises the interburst-interval (IBI) signal.Based on visual examination, it could be observed that the bursting firing behaviour in vagal signals was consistently aligned with respiration patterns of inspiration-expiration.As a result, the breathing rate (BR) was extracted from the bursting activity by expressing the IBI signal in bursts per minute (see supp.figure 3 for an overview of the methodology).The raw signal underwent three decomposition strategies to get the inter-burst-interval signal from the power spectrum, cardiac peaks from low-frequency field analysis and continuous wavelet analysis, and neural peaks from high-frequency analysis and SO-CFAR threshold.This analysis allows us to get (d) the breathing rate, (e) the heart rate HR, and (f) the metrics evolutions from the clusters.(g) Glucose levels were collected every 5 min.For each observation (i.e.recordings), (h) the multivariable correlation matrix was computed using these signals as inputs.Percentages were then computed within cohorts and observation types (glucose-increasing or decreasing recordings) to summarise the population results.LPF: low-pass filter, BPF: band-pass filter, HR: heart rate.

Neural component 2.3.4.1. Pre-processing
The raw signals acquired at 30 KHz were first downsampled to 10 KHz to preserve the original signal while optimising the subsequent processing.Then signals were filtered using a 9th order Butterworth, this time using a bandpass filter with cut-off frequencies at 400 Hz and 5 KHz to extract the high frequencies of interest.A second-order infinite impulse response (IIR) notch digital filter was also applied at odd harmonics of 50 Hz between 350 Hz and 5 KHz to reduce line noise from electronic current sources.The activation of the vagus nerve in response to blood glucose levels was analysed by dividing each 120 min recording into regions of at least 20 min duration containing a clear increase or a clear decrease in blood glucose, with each region constituting one observation.By selecting periods of significant change of at least 20 min, we could observe clearer patterns and relationships, which helped to better analyse the physiological and neural regulatory mechanisms involved.Either none, one or more than one of such observations could be recognised throughout an experimental protocol.Each of the observations underwent the steps described in the following sections (further detailed in supp.figures 4 and 5).

Thresholding and peak detection
The amplitude of the neural recordings is typically not stable over time due to changes in the impedance of the electrode-tissue interface.In these recordings, we also found an amplitude-modulation of the signal due to breathing.For this reason, a thresholding method was implemented based on the smallest of constant false-alarm rate (SO CFAR) filter (Scharf 1991), which has been successfully implemented in previous studies to accommodate for the respiratory modulation on vagus nerve recordings (Zanos et al 2018).CFAR filters use a sliding window to continuously estimate the background mean power of a signal so that a threshold that maintains a CFAR on average can be applied on a per-window basis.This allows the threshold to adapt to changes in the background mean power.The parameters of the threshold that were heuristically chosen and used in all recordings were the window duration (150 ms on each side), the guard cell duration (1 ms on each side), and the threshold level (3 SDs from the mean).Windows of 8 ms centred around peaks surpassing the threshold were extracted.Peaks with a width of less than 0.2 ms (between zero-crossing points), or an amplitude of more than five times SDs from the mean were discarded.In addition, peaks detected within an 8 ms window centred around a cardiac peak were also excluded.

Dimensionality reduction
The high-dimensional space corresponding to the waveforms extracted after cardiac waveform elimination was reduced using linear and nonlinear dimensionality reduction methods.In particular, uniform manifold approximation and projection (UMAP) was applied.UMAP is a non-linear method that outperforms the linear principal component analysis (PCA) in datasets where non-linear dependencies are remarkable (Nanga et al 2021, Xiang et al 2021).Dimensionality reduction was necessary as modelling with high-dimensional raw signal data may lead to overfitting and reduce the clustering performance.The dimensionality of the original space, 40 samples corresponding to the length of the detected compound action potential (CAP) waveform, was reduced to three dimensions.This number of dimensions was heuristically chosen and guided by previous literature (Zanos et al 2018).Data was first standardised (zero mean and unit variance) to ensure all features were at the same scale and preserve the structure of the original data.

Clustering
The main hypothesis for sorting CAP waveforms into different groups (clusters) is that different CAP waveform shapes and amplitudes represent the activity of different groups of nerve fibres.Once dimensionality reduction was performed, two unsupervised clustering algorithms were implemented.The K-means algorithm was generally used to keep consistency among trials on the number of identified clusters because this algorithm requires the number of clusters to be specified and was heuristically set to four after preliminary analysis to keep a trade-off between clusters being significantly different in all recordings and reducing the computational power.
The algorithm clusters data by trying to separate samples into n groups of equal variance, minimising a criterion known as the inertia or within-cluster sumof-squares.On a subset of shorter recordings (10-15 min duration), the density-based spatial clustering of applications with noise algorithm was used to automatically identify the maximum number of distinct clusters.The algorithm finds core samples of high density and expands clusters from them.The two parameters of the method were chosen heuristically after a preliminary analysis with reference to a previously published methodology for the sake of consistency and comparison (Zanos et al 2018): the maximum distance between two samples for one to be considered as in the neighbourhood of the other was set to 0.5, and the minimum number of points required to form a cluster was set to 30.

Neural decoding: metrics computation
Each cluster of CAP waveforms was inspected visually by the average waveform.The distribution of inter-spike-intervals (ISIs) was computed to reveal the nature, noise or neural, of the recorded spikes.To obtain the metrics for the detected clusters, a rolling window of 10 s was used to compute the spike rate and the mean amplitude of the peaks within each window for each cluster (figure 1(f)).These metrics were also extracted for all the detected peaks to get information about the overall activation of the nerve.

Correlation
For each observation in each animal (i.e.glucoseincreasing or glucose-decreasing period), pairwise correlation using the Pearson standard correlation coefficient (ρ) was conducted among all neurallyextracted metrics (i.e.spike rate and the mean amplitude of the peaks) and physiological parameters (BR extracted from IBI signal, HR, and blood glucose concentration).The glucose measurements occurred at a lower frequency (every 5 min) compared to the neural-extracted metrics (10 KHz sampling frequency).To combine both temporal scales and determine the correlation coefficient computations, the average of the neural metrics and non-glucose physiological parameters were obtained from a 5 min window prior to each glucose measurement.The 5 min window was chosen to avoid interpolationrelated issues on the glucose signal and obtain a representative trend of the other signals to enable a more comprehensive analysis of linear relationships.To provide more informative results, the resulting correlation matrix was plotted for each observation showing the coefficients on the top right side, scatterplots showing the linear regression between the two compared variables on the bottom-left side, and the distributions for each individual variable on the diagonal (figure 1(h)).Observations and animals were defined as responders when a strong positive or negative correlation was found in at least one cluster, or a weak correlation was present in the majority of the clusters.Correlation was defined as strong when |ρ| ⩾ 0.7, weak when 0.5 ⩽ |ρ| ⩽ 0.7, and no correlation |ρ| ⩽ 0.5.The percentage of responder (strong or weak) and non-responder observations (or animals) in each experimental condition (intact, afferent, or efferent) and observation type (glucose-increasing and decreasing) was calculated by dividing the number of observations with strong, weak and no correlation coefficients by the total number of observations for each observation type in that condition, expressed as a percentage.All the above-mentioned analysis was conducted using python and its relevant libraries.

Results
In this work, we used gold hook microwires to record high signal-to-noise neurograms from the surface of the cervical vagus nerves of anaesthetised Sprague-Dawley rats (n = 11 in intact cohort, n = 6 in each afferent and efferent recording cohort) (Güemes Gonzalez et al 2023).The noise floor of all recorded signals was 6.4 ± 1.3 uV.The average impedance for all recordings was 65.9 ± 36.9 kOhm.No correlation was found between the noise floor and impedance values in any of the groups (intact group: ρ = 0.34, ligated group: ρ = −0.29).In addition, we observed similar impedance values in both responders and non-responders, indicating that the observed effects are not solely reliant on impedance variations.More details on the background noise and impedance values for each cohort can be found in supp.figure 6.

Decoding physiological signals from vagus nerve recordings
The raw recording contained a variety of signals, with some arising from non-neural sources like cardiac and skeletomuscular activity, which was identified and eliminated as described in sections 3.1 and 3.2.
In addition to the cardiac and muscular contribution, in these recordings, we found a clear amplitudemodulation of the signal on the temporal domain that corresponded to respiration through visual inspection.This relation to breathing patterns was even more remarkable in the frequency domain, where bursts of high-frequency spikes could be observed (figure 1(c)).To validate this observation, a small pilot study (n = 3) was conducted under inhaled anaesthesia to control the BR by adjusting the concentration of isoflurane, which ranged from 5% to 1.5%.As expected, the periodicity of the bursts of amplitude modulation was always aligned with the BR of the animal (figure 2).Different strategies were implemented to understand if the increased amplitude found during these bursts was due to micromotions of the vagus nerve during breathing, which could account for changes in the impedance of the electrode-nerve interface, or if these occur as a result of neural activity.The breathing modulation was absent in the time-domain recordings obtained from the second electrode placed on muscles in the proximity of the vagus nerve, and the contribution in the frequency domain was clearly diminished, with bursts not exceeding 700 Hz.Secondly, the application of lidocaine eliminated or significantly attenuated the amplitude of these bursts to up to 30% of the original activation, both in the temporal and the frequency domains, where the response became similar to the recordings from the muscle (figure 2(b)).These reductions were more significant in the intact and afferent recordings (maximum frequency contribution of the pre-lido breathing-bursts exceeded 3000 Hz) than in the efferent recordings (maximum burst frequencies of 1600-3000 Hz).The BR was therefore extracted from the IBI signal calculated using the sum of the power of frequencies over 500 Hz.
The contribution of the breathing modulation to the neural recording was reduced using an adaptive thresholding method based on the SO CFAR filter that was successfully implemented in previous work to accommodate for the respiratory modulation on vagus nerve recordings (Zanos et al 2018).By following this methodology, we were able to identify and extract neural spikes that occur in between and during the respiratory modulations.The identified spikes presented a variety of waveforms, the majority of them being asymmetrical (supp.figure 7(a)).The exact waveform of the spikes detected from neural recordings, i.e. shape and amplitude, is representative of different characteristics of the nerve activation, such as the number of fibres discharging simultaneously, and the size and location of the fibres relative to the recording electrode placed at the surface (Parker et al 2018, Chang et al 2019).As a result, these spikes can be considered compound neural action potentials (CNAPs) that arise from the contribution of multiple single units (Zanos et al 2018, Chang et al 2019, Koh 2019, Masi et al 2019).All identified neural spikes were clustered into groups under the hypothesis that different clusters correspond to different activation patterns.UMAP allowed to reduce the high-dimensional space corresponding to the extracted waveforms to only three dimensions, which increased the efficiency of the clustering algorithms by removing redundancy in the data.
To verify that the detected spikes were of neural nature, a lidocaine bath at the exposed site used at the end of the experiments confirmed that the response was abolished or significantly attenuated.Moreover, neural spikes show an ISI of 9 ms to 11 ms as previously reported (Jiman et al 2020), whereas noise spikes appear with a smaller ISI of 4 ms, coinciding with harmonics of the 50 Hz electric noise (supp.figure 7(b)).The distribution of ISI times for the detected spikes was therefore used as an additional indicator of the nature of the recordings.

Glucose decrease, but not increase, is consistently conveyed by the intact vagus nerve
In all the animals anaesthetised with urethane we observed an increase in glucose levels during the first 30 min due to an increase in sympathetic outflow driven by this anaesthetic (Sánchez-Pozo et al 1988) and the surgical stress, which was followed by a decrease in glycaemia.Urethane-induced hyperglycaemia was observed again 60 min post-injection, generally coinciding with the 1 ml IP injection of saline at minute 30 (figure 3).A 1 ml IP dose of insulin injected 20-50 min after the saline injection   a progressive continuous reduction glucose levels from 5 to 20 min after injection independently of the vagotomy.In a total of four animals with ligated nerves (two afferent and two efferent), an additional IP injection of glucose (2 g kg −1 ) was done to evaluate differences with urethane-induced hyperglycaemia and the impact of injection timing (figure 1(a)).
The activity of the intact vagus nerve in response to changes in glucose levels was assessed by identifying and clustering CNAP from 30 min periods around each IP injection (10 min before and 20 min after injection).We observed an increase in the overall firing rate and in the firing rate of specific clusters after insulin injection corresponding to a decrease in glucose levels, which was accompanied by a decrease in the amplitude of the detected spikes (figure 3(b)).No consistent changes in the firing rate were observed on increasing blood glucose concentrations (figure 3(a)).
The activation of the intact vagus nerve in response to blood glucose levels was further analysed by dividing each recording into regions of at least 20 min duration containing a clear increase or clear decrease in blood glucose, with each region constituting one observation.More than one of such regions could be recognised throughout an experimental protocol (e.g.there are a total of seven glucose-increasing observations in six animals showing increases in glucose levels).An illustrative example of the obtained pair-wise correlation matrix used to confirm and quantify the degree of the linear relationship between the neural metrics, and the physiological events (glucose, ECG and breathing) can be observed in figure 4.An overview of the correlation coefficients for all animals and observation types is available in supplementary tables 1-4.The percentage of strong, weak and decorrelated pairs within each observation type (glucose-increasing and decreasing) and experimental condition (intact, afferent, or efferent) can be found in table 3.
In the intact cohort, 91% of the animals presented a neural change due to glucose decreases, in the form of increased firing rate (100% of the time) and amplitude (90% of the time, table 1).In the majority of the intact nerve observations corresponding to a decrease in glucose levels, we found a strong increase (ρ ⩽ −0.7) in the overall firing rate and in most of the responsive clusters with glucose.Similarly, in most of such observations, we found an increase in the amplitude of certain clusters corresponding to a reduction in glucose levels (ρ ⩽ −0.5), but no conclusive changes were observed in the overall nerve signal amplitude.On the other hand, when glucose concentrations increased, we found no convincing correlations (|ρ| ⩽ 0.5) of the firing rate or the amplitude with glucose, with a small percentage of the observations presenting a negative correlation with increased glucose.At the animal level, however, 60% of the responsive animals showed a decrease in the firing rate (table 2).

Vagus nerve firing rate associated with decreases and increases in blood glucose is conveyed by both afferent and efferent pathways
When recording afferent sensory signals from the distal terminal, we found a negative correlation of glucose with the overall firing rate and most clusters in the majority of glucose-decreasing observations (ρ ⩽ −0.5) (figure 5(a)).In recordings with increasing glucose concentrations, we found no meaningful correlations (|ρ| ⩽ 0.5) of the firing rate with Table 1.Percentage of animals for each independent experimental protocol (intact, afferent and efferent) responding to glucose-decreasing changes.Responders are considered when there was a strong positive or negative correlation in at least one cluster, or a weak correlation was present in the majority of the clusters.Note that only in five animals with efferent nerve ligation (out of 6) there was at least one clear glucose-decrease observation.glucose in most of the observations, but responsive clusters showed an increased firing activity in half of the observations (ρ ⩾ 0.5).Glucose was negatively correlated (ρ ⩽ −0.5) with the amplitude of the overall nerve activation in the majority of the glucose-decreasing observations.Interestingly, the afferent responsive clusters reduced their amplitude with glucose in slightly over half of the glucosedecreasing observations and increased their amplitude in slightly below half of glucose-increasing observations (ρ ⩾ 0.5).
When recording motor efferent signals from the proximal terminal we found no correlation between the overall activation and glucose in most of the glucose-decreasing observations (|ρ| ⩽ 0.5), but we found at least one efferent responsive cluster increasing the firing rate corresponding to decreases in glucose in a similar number of the observations (ρ ⩽ −0.5, figure 5(b)) and in a larger number of observations where glucose was rising (ρ ⩾ 0.5).In the majority of the cases, these changes in the firing rate of responsive clusters were also accompanied by a decrease in their amplitude in glucose-decreasing observations and an increase in their amplitude in glucose-increasing observations (ρ ⩾ 0.5 in either case).An illustrative example of the temporal evolution of metrics in the afferent and efferent recordings is depicted in figure 5, and the corresponding correlation matrix in supp.figure 9.

Amplitude modulations are mostly driven by respiratory and cardiac signals
We further evaluated the possibility that the changes in firing rates and amplitude evolutions were due to changes in breathing by measuring the correlation of all variables with the BR signal obtained from the IBIs signal (table 3).We generally found no conclusive correlations between the BR and firing rate in intact and efferent observations independently of the glucose trend.However, we found BR to be weakly correlated with the firing rate in most glucose-increasing and more than half of glucose-decreasing afferent recordings (ρ ⩾ 0.5).The contribution of respiratory signals to amplitude modulation was more remarkable accounting for most of the glucose-increasing observations in the intact nerve (ρ ⩾ 0.7), and around 70% of afferent recordings (ρ ⩾ 0.5), and efferent recordings (ρ ⩽ −0.5).On glucose-decreasing observations, results in the intact and efferent nerve were not consistent, but were positively correlated with glucose decreasing in a large percentage of afferent recordings (ρ ⩾ 0.5).No meaningful correlation was found between BR and the other physiological variables, HR and glucose.
Regarding HR, we found it to be positively correlated with the overall nerve signal amplitude (ρ ⩾ 0.5) in most of the glucose-increasing and decreasing observations in the intact cohort.In glucose-decreasing observations, HR was majoritarily found to be weakly positively correlated (ρ ⩾ 0.5) to the amplitude of the nerve activation in the efferent recordings, with no correlation found in the afferent recordings.In glucose-increasing observations, no correlation was found between HR and the amplitude of the nerve activation in most of the efferent observations, whereas in the afferent recordings a positive correlation was frequently observed with amplitude (ρ ⩾ 0.5).Similarly to the respiratory signalling modulation, no significant correlation of HR was found with firing rate in any of the cohorts, independently of the glucose trend.
In terms of the relationship between HR and glucose, we frequently found a weak positive correlation between glucose and HR in recordings from the intact nerve when glucose was decreasing, which arises from the consistent decrease in the trend that the HR exhibited in all the intact nerve experimental protocols (supp.figure 8).Glucose was also weakly negatively correlated to HR in some of the glucose-increasing observations in the intact nerve (ρ ⩽ −0.5), showing no correlation in the remaining observations.This further confirms the progressive decrease in HR over time.No correlations were found in the majority of the afferent recordings, when glucose was increasing, and a few of the efferent recordings showed a negative correlation.When glucose was decreasing, HR weakly correlated to glucose in less than half of efferent and afferent observations (ρ ⩾ 0.5).
The positive correlation obtained between HR and glucose required further evaluation to determine if the observed changes were related to glucose or HR.A similar methodology was followed dividing the recordings this time into observations based on increasing and decreasing HR trends.This approach reduces the impact of the observations selected previously based on glucose trends.In observations from afferent recordings where the HR trend was clearly increasing, no correlation was obtained with the amplitude and the firing rate of the vagus nerve in all the observations, nor with glucose in the majority of them.Similarly, in HR-increasing observations from efferent recordings, no correlation was found with the firing rate (all the observations) or glucose (half of the observations), but a weak positive correlation with amplitude was frequently observed.No meaningful correlation was found with any of the metrics in efferent HR-decreasing recordings.There was a strong negative correlation with firing rate accompanied by a strong positive correlation with glucose in most of the afferent recordings, as observed previously.See supp.table 5 for a detailed overview of the correlation coefficients.

Discussion
Understanding the transmission of metabolic signals by the vagus nerve in isolation is a significant Table 3. Overview of population percentages of correlation coefficients for all pair-wise analysis conducted in each experimental cohort and observation type.Responsive observations show at least one cluster showing strong correlation, or a weak correlation present in the majority of the clusters.Strong correlation: |ρ| ⩾ 0.7, weak correlation: 0.5 ⩽ |ρ| ⩽ 0.7, no correlation: |ρ| ⩽ 0.5.In bold convincing outcomes when more than 60% of the observations met the same outcome.Note the number of recordings for each observation type (glucose trend) differs from the number of animals as there could be more than one or no observations for a certain type.challenge because it is involved in many physiological processes, including digestion, HR, and immune response.This paper presents for the first time a multivariable physiological model of vagus nerve transmission (table 4) and a framework for studying signalling mechanisms in peripheral autonomic nerves.

Intact
A spike-sorting methodology is proposed for recordings from the surface of the nerve, and a comprehensive multivariable correlation analysis is conducted.The metrics derived from the recorded vagal nerve activity, firing rate and amplitude, showed intriguing patterns linked to BR, HR, and blood glucose concentrations.The evolution of the amplitude of the detected spikes over time can provide insightful information about the contribution of the population of fibres in the nerve.However, careful consideration is needed when analysing this metric to discard non-neural changes (e.g.motion artefacts and sharp changes in impedance).
All the recordings presented a marked breathingmodulation both in the temporal and frequency domains.In this work, our strategies based on nerveblocking using lidocaine and recording from nonneural regions revealed the neural nature of these signals, but do not fully discard per se a potential amplification effect arising from transitory changes in impedance due to breathing motion.Previous work has confirmed that respiratory information is conveyed by specific fibres of the vagus nerve, as demonstrated using surface (Zanos et al 2018, Jiman et al 2020, Marmerstein et al 2021, Vallone et al 2021) and intraneural electrodes, where modulation was present only in few contacts (Jiman et al 2020).In this work, we present a novel strategy based on IBI that allows to easily extract the temporal evolution of the breathing-modulated bursts and the BR over time and differs from the root mean square or the power of the recorded signal as previously used (Sevcencu et al 2018).Our findings suggest that the afferent signalling related to breathing drives the amplitude and firing rate modulation, while efferent signalling varies depending on glucose trend.Sevcencu et al (2018) demonstrated a strong correlation between the temporal profile of respiratory signal derived from the vagus nerve and the respiratory signal obtained using a pressure transducer.This similarity is due to the encoding of respiratory information by large myelinated afferent fibres from pulmonary stretch receptors in the lungs, with maximum recruitment achieved at the end of inspiration (Ottaviani et al 2020, Patros et al 2022).We similarly demonstrated that the amplitude of the breathing modulation of the afferent and intact recordings is larger than in efferent recordings, whose activity could originate from the laryngeal fibres that fire motor action potentials during lung inflation (Green andNeil 1955, Sevcencu et al 2018).Modulations due to changes in HR in the intact and efferent recordings show parasympathetic control, with the motor vagus nerve activating to lower HR in tachycardia (increase in efferent amplitude) and increasing sensory activity following bradycardia (increase in afferent firing rate) (Capilupi et al 2020, Patros et al 2022).
In this study, we have primarily focused on the specific contribution of the vagus nerve to the transmission of glucose-related signals.Recording from the cervical vagus nerve account for signalling corresponding to multiple internal organs.Therefore, in this discussion, we analyse how the individual vagal branches of different internal organs respond to changes in glucose to dive into the physiological origin of our results.We found a progressive decrease in the firing rate of the intact vagus nerve and afferent fibres in 60% and 40% of the responsive animals respectively as glucose levels were rising due to the effect of urethane.This has been supported by previous studies in isolated and perfused livers, and in situ animal models, where glucose-inhibited vagal afferents within the hepatic branch of the vagus nerve reduced their firing activity with increases in the portal vein glucose concentration (Niijima 1982, Niijima 1984, Yamatani et al 1998).This mimics the situation that occurs after feeding when the glucose level is higher in the portal vein, than in the arterial circulation, creating a negative arterial-portal glucose gradient (Jackson et al 1997).Our recordings from the cervical nerve, however, show a larger variability in the sensory responses to hyperglycaemia than when recording directly from the branches of interest, which aligns with previous studies recording from the surface or within the cervical vagus nerve following IP glucose injection (Masi et al 2019, Jiman et al 2020).One reason may be due to the differential contribution of afferent branches from different organs.For example, studies on gastric motility have shown that peripheral hyperglycaemia stimulates vagal afferent pathways from the gastroduodenal mucosa to the solitary tract (NTS) and the dorsal motor nucleus of the vagus to mediate gastric relaxation (Zhou et al 2008), which is contrary to the inhibition of activity observed the hepatic branch (Niijima 1982, Niijima 1984, Yamatani et al 1998).Additionally, the use of urethane may have reduced the negative arterial-portal gradient created by normal feeding, intra-portal, or intravenous (IV) administration of glucose, therefore reducing the response of afferent fibres the response to hyperglycaemia.A larger animal cohort will be needed to provide more conclusive results in this regard.On the other hand, when recording from the efferent branch, we found clusters with an increased firing rate and increased amplitude of the detected compound action potential in over 75% of the recordings, which relates to a higher number of fibres being simultaneously active.This is also in line with studies where elevated portal vein glucose after IV administration increased the efferent firing rate of the hepatic branch inhibiting glucose production (Niijima 1985), and of the motor pancreatic branch of the vagus to stimulate insulin secretion (Yamatani et al 1998, Donovan andWatts 2014).The responses observed due to activation of vagal efferent branches have also been confirmed using vagus nerve stimulation and vagal blocking (vagotomy or use of tropine) in several species including rats, pigs, dogs and humans (Blat andMalbert 2001, Meyers et al 2016).
In insulin-induced hypoglycaemia, we found a progressive increase in firing rate in the intact vagus nerve (91% of the animals in the intact cohort), which aligns with previous work using cuff electrodes placed on the vagus nerve in mice (83.3% of experiments injecting insulin) (Masi et al 2019), and slightly exceeds the results reported in rats using intraneural carbon-fibre microelectrode array electrodes that depends on the specific fibres that are being interrogated (66.7% of the experiments) (Jiman et al 2020).We similarly found an increase in the firing rate of afferent clusters in 72% of the recordings (100% of the afferent responsive animals), which was accompanied by an increase in the amplitude of the detected cluster.Glucose-sensitive vagal hepatic afferents appear to have no role in mediating hypoglycaemic detection at the portal-mesenteric vein based on experiments blocking nerve transmission through vagotomy or vagal cooling (Niijima 1982, 1984, 1985, Donovan and Watts 2014).However, sensing of low glucose concentration in the portal vein following hyperglycaemia by glucose-inhibited hepatic afferents in the liver could have facilitated the firing activity of the hepatic vagus nerve.Moreover, injections of insulin in the isolated and perfused liver preparations (Niijima 1983), and from intraneural recordings in the cervical vagus nerve, have reported increased afferent activity (Jiman et al 2020), which may resemble the paradoxical signalling of hunger produced by endogenous insulin after meal intake (Niijima 1969, Rezek et al 1975, Rezek 1976).We found that the increase in the intact nerve activity led by glucose-decreasing trends is not solely driven by afferent activity, with motor clusters also showing an increased firing rate in 63% of the recordings.In this line, Masi et al found that the ability of a regression model to recover blood glucose from vagus nerve recordings was significantly reduced in transient receptor potential vanilloid subfamily V member 1 (TRPV1) cell depleted experiments compared to the wild-type animals (average regression error increased from 5.1 mg dl −1 to 18.3 mg dl −1 ), but the algorithm was still capable to recover the decreasing profile of blood glucose levels suggesting other signalling mechanisms being involved.TRPV1 receptors largely contribute to regulating appetite hormones and increasing the activity of gastrointestinal vagal afferents after meal intake (Suri and Szallasi 2008, Zsombok 2013, Derbenev and Zsombok 2016).An increase in motor firing rate during insulin-induced hypoglycaemia is largely supported by previous evidence of a central regulation to coordinate the function of peripheral organs (Wu et al 2004, Pocai et al 2005a, 2005b, DiCostanzo et al 2006), as opposed to the peripheral regulation observed in hyperglycaemia.Our findings of increased motor activity can relate to any of the signalling mechanisms triggered by insulininduced hypoglycaemia, such as stimulation of pancreatic efferent nerve activity by orexin neurons located in the lateral hypothalamic area (Wu et al 2004) or by glucose transporter 2 (GLUT2) neurons in the NTS (Lamy et al 2014) that drive the counterregulatory response to hypoglycaemia by increasing glucagon secretion, or increased efferent vagal drive to the stomach to promote gastric emptying (Hjelland et al 2005, Zhou et al 2008).Increased insulin sensing directly in the arcuate also reduces hepatic glucose production via activation of the motor vagus nerve (Pocai et al 2005a(Pocai et al , 2005b)).It is worth mentioning that the timing and dynamics for observing changes in nerve activity in response to variations in blood glucose level depend on the animal species and the specific dose and route of insulin or glucose administration.For example, Niijima (1983) showed changes in afferent activity after 5 min of intra-portal venous administration of 20 mU of insulin in the guinea pig, and a few years later (Niijima 1985) he reported an immediate decrease in the efferent activity of the nerve after 4 u IV injection of insulin at the jugular vein in rats.In our study, the neural activity response was observed promptly upon detecting a decrease in glucose levels, typically with a delay of around 7.5 min following IP injections.These findings align with similar observations reported by Jiman et al (2020) also in rat models.Conversely, in mice, Masi et al (2019) reported an immediate change in glucose levels following IP injection, accompanied by immediate responses in vagus nerve activity, likely due to the smaller size of this particular model.The contribution of other physiological parameters derived from the neural recordings (HR and breathing) to the changes in the firing rate was not meaningful from our recordings.The positive correlation of glucose and HR during glucosedecrease recordings evidences the consistent decrease in HR over time because of urethane, which has been previously reported to lower significantly resting HR and systolic blood pressure as compared to unanaesthetised animals (Kubota et al 1756, Maggi andMeli 1986).
An intrinsic limitation of any terminal work is the use of anaesthetics (Jiman et al 2020).In this work, urethane was chosen as the preferred anaesthetic because it does not depress neural activity as much as some common alternatives such as isoflurane (Maggi andMeli 1986, Hemmings et al 2005).However, it has been shown to have a large impact on glucose metabolism by increasing sympathetic activity and affecting the action of exogenous insulin, resulting in a physiological scenario that is far from the observed in awake animals (Wang et al 2016).Chronic experiments providing a physiological evaluation of the vagal signalling throughout the full digestive process are strongly encouraged as they will enable a better understanding of the vagus transmission and confirm our conclusions.In the recovery setup, the translation of microwires is hindered by a number of limitations, such as the need for good encapsulation and stability in chronic implantation.Future work is therefore needed to evaluate the use of flexible and conformable devices as a solution for long-term studies.Another limitation of this work is the lack of ability to continuously record HR and breathing using independent monitoring systems, which will fully isolate these signals and facilitate a detailed analysis of the time association between the physiological change and the observed neural response.We, however, have successfully identified this information from the neural signals and consistently verified their relationship with the physiological processes of interest during the experiments.Additional studies are also encouraged to keep validating these findings on a larger population and contribute to a more comprehensive understanding of the phenomena being investigated.Finally, in this work have derived our conclusions from a correlation analysis, which fails to identify non-linear relationships (e.g.sudden rather than progressive changes) and does not imply causation.Other statistical methods that capture non-linear and casual relationships, such as structural equation modelling, are encouraged in the future for modelling the complex relationships observed between multiple variables.As a result, the conclusions presented in this study should be evaluated bearing in mind the experimental and analytical methodology that has been described.
Closed-loop bioelectronic medicine shows promise for diagnosing and treating chronic diseases like diabetes.Advances in neural decoding algorithms and data analysis techniques allow to identify information about the state of nerves and related organs (Vitale and Litt 2018, Zanos 2019).The work presented in this paper progresses towards leveraging the metabolic information decoded from the vagus nerve to improve glucose control in people with diabetes.Fast detection of neural signals related to changes in blood glucose is crucial to preventing risky hypoglycaemic or hyperglycaemic events.Purely metabolic strategies, like the current artificial pancreas, disregard the important contribution of the nervous system in healthy glucose regulation leading to suboptimal performance (Güemes and Georgiou 2018, Güemes Gonzalez et al 2020).Future work will evaluate the benefits of integrating the information from neural recordings and biomarkers to obtain faster and more accurate information about the metabolic state, which can then be used to determine optimal insulin doses and deliver precise electrical stimulation to restore metabolic balance.Such a neuro-metabolic closed-loop platform has the potential to enhance patient autonomy, minimise delays in insulin action, and reduce the risk of misestimating carbohydrate content during meal intake.

Conclusion
In this work, we have presented a framework to decode the motor and sensory vagal signals of glucose, HR and breathing using signal processing and descriptive data mining techniques, including filtering, wavelet decomposition and correlation analysis.We have also introduced a novel strategy to extract the BR over time based on IBI.This analysis has contributed to the development of a multivariate physiological model of vagus nerve transmission by providing important insights into the complex interactions between the central nervous system, the autonomic nervous system, and the metabolic system.This model advances our understanding of how the vagus nerve encodes information to control glucose metabolism and homeostasis.Ultimately, this work establishes a promising new pathway for the development of closed-loop neuro-metabolic systems that can exploit neural control for improved diabetes diagnosis and treatment.

Ethical statement
All animal work was carried out in accordance with the UK Animals (Scientific Procedures) Act, 1986.

Figure 1 .
Figure 1.Overview of methodology.(a) Timeline of experimental protocol (in blue animals undergoing nerve ligation).(b) Raw signals were recorded from the cervical vagus nerve.(c)The raw signal underwent three decomposition strategies to get the inter-burst-interval signal from the power spectrum, cardiac peaks from low-frequency field analysis and continuous wavelet analysis, and neural peaks from high-frequency analysis and SO-CFAR threshold.This analysis allows us to get (d) the breathing rate, (e) the heart rate HR, and (f) the metrics evolutions from the clusters.(g) Glucose levels were collected every 5 min.For each observation (i.e.recordings), (h) the multivariable correlation matrix was computed using these signals as inputs.Percentages were then computed within cohorts and observation types (glucose-increasing or decreasing recordings) to summarise the population results.LPF: low-pass filter, BPF: band-pass filter, HR: heart rate.

Figure 2 .
Figure 2. Understanding the neural nature of breathing modulation.(a) Pilot study using isoflurane to change the breathing rate and visualise the relationship with the neural bursts detected in the spectrogram, where the temporal distance between bursts increases with deeper isoflurane concentrations.(b) Reduction of burst activity in the spectrogram when recording from locations in nearby muscle and when applying lidocaine to the nerve, compared to the high-frequency bursts of activity recorded when placed around the nerve.(c) Lidocaine was used at the end of each experiment to evaluate the neural origin of the recorded signals.A clear drop in the firing rate and amplitude of the signal was observed after the application (grey line).

Figure 3 .
Figure 3. Illustrative example of one animal in the intact nerve cohort showing the temporal evolution of breathing rate, overall spike rate and mean amplitude of the detected clusters, and the firing rate of each of the four clusters corresponding to increases and decreases of glucose level.Superimposed smoothed signals obtained from a moving average using 5 min windows are shown to facilitate the visualisation of trends.Inserts in the left-and right-hand sides of the image depict the average waveform of each identified cluster and the histogram corresponding to the inter-spike-interval (ISI) metric.(a) Glucose-increasing recording: the overall firing rate fluctuates over time, whereas the amplitude of detected spikes clearly increases, and the firing rate of clusters shows a range of responses (the blue cluster decreases, the green increases, the orange slightly increases, and the purple fluctuates).The 'disconnection' label refers to a transient mechanical artefact.(b) Glucose-decreasing recording: the overall firing rate increases as glucose starts to drop, with the amplitude of detected spikes decreasing over time.The firing rate of most of the identified clusters (blue, orange and green) shows an increase over time as glucose drops.

Figure 4 .
Figure 4. Correlation matrix corresponding to recordings from the same animal (intact nerve cohort) in figure 3. Correlation matrix corresponding to the glucose-increasing observation (signals depicted in figure 3(a)).(b) Correlation matrix corresponding to the glucose-decreasing observation (signals depicted in figure 3(b)).

Figure 5 .
Figure 5. Illustrative example of one animal in the (a) afferent cohort and (b) another animal in the efferent cohort showing the temporal evolution of breathing rate, overall spike rate and mean amplitude of the detected clusters, and the firing rate of each of the four clusters corresponding decreases of glucose levels.Superimposed smoothed signals obtained from a moving average using 5 min windows are shown to facilitate the visualisation of trends.Horizontal lines in the overall spike rate depict maximum and minimum values in smoothed signals for further visualisation of the increasing trend.Inserts in the left-and right-hand sides of the image depict the average waveform of each identified cluster and the histogram corresponding to the inter-spike-interval (ISI) metric.The sharp change in signals in panel (b) around minute 65 corresponds to a sudden and strong contraction of neck muscles, which despite being detected and eliminated to the greatest extend, still interfered with this recording.The purple cluster in the efferent recording (bottom right) illustrates a clear non-neural signal in both the waveform and the ISI distribution and was not considered in the correlation analysis.The corresponding correlation matrix can be found in supp.figure 9.

Table 2 .
Percentage of animals for each independent experimental protocol (intact, afferent and efferent) responding to glucose-increasing changes.Responders are considered when there was a strong positive or negative correlation in least one cluster, or a weak correlation was present in the majority of the clusters.Note that only in six animals with the intact vagus nerve (out of 11) there was at least one glucose-increase observation.

Table 4 .
Multivariate physiological model of the vagus nerve modulation of firing rate and amplitude due to metabolic, cardiac and respiratory changes.