Browse the corpus

A probabilistic model of behavioural emergence from general anaesthesia in mice. BACKGROUND: Time to emergence from general anaesthesia is highly variable between individuals. This variability has been attributed to individual differences in anaesthetic sensitivity. However, this hypothesis has not been verified experimentally. We explicitly test this hypothesis by quantifying emergence from anaesthesia repeatedly in the same individuals over time. METHODS: Genetically identical adult (12-24 weeks old) male (n=40) and female (n=20) C57BL/6J mice were exposed to 2 h of isoflurane (0.90 vol%) on 10 separate occasions. Time to emergence was measured using the return of the righting reflex. Predictions of the standard effect-site pharmacokinetic-pharmacodynamic (PK-PD) model and neuronal dynamics model of stochastic fluctuations between the awake and anaesthetised states were fit to observed emergence times. Repeated steady-state assessments of the righting reflex obtained during the last 2 h of a 4-h exposure to 0.3, 0.4, 0.6, or 0.7 vol% isoflurane (n=20 per concentration) were used to determine individual probabilities of losing the righting reflex, which was defined as an individual's anaesthetic sensitivity. RESULTS: Emergence times varied by at least two orders of magnitude after identical anaesthetic exposure. We did not find consistent inter-individual differences in emergence times. Instead, we found that variability in emergence times across trials in each individual was as large as that between two different individuals. Emergence times were not correlated across time. Consistent with previous work, we identified large individual differences in anaesthetic sensitivity which persisted on a time scale of at least 1 week. A standard PK-PD model failed to reproduce inter-trial variability. In contrast, the neuronal dynamics model reproduced both population- and individual-level variability in emergence times. CONCLUSIONS: Stochastic state switching contributes to inherent variability in emergence from general anaesthesia. Delayed emergence occurred in a small proportion of anaesthetic exposures in a genetically homogeneous population. The neuronal dynamics model predicts that anaesthetic emergence times will be probabilistically long, which might explain delayed emergence observed in clinical settings.

Studies were approved by the Institutional Animal Care and Use Committee at the University of Pennsylvania, and were conducted in accordance with the US National Institutes of Health guidelines. Inbred male and female wild-type C57BL/6J mice (inbred generation F226 or higher, >99.99% genetic identity28; Jackson Laboratory, Strain 000664, Bar Harbor, ME, USA) aged 12–24 weeks were used to study anaesthetic emergence (n=40 males, 20 females) and anaesthetic sensitivity (n=60 males).

titutes of Health guidelines. Inbred male and female wild-type C57BL/6J mice (inbred generation F226 or higher, >99.99% genetic identity28; Jackson Laboratory, Strain 000664, Bar Harbor, ME, USA) aged 12–24 weeks were used to study anaesthetic emergence (n=40 males, 20 females) and anaesthetic sensitivity (n=60 males). Before anaesthetic exposures, mice were habituated to gas-tight, 37°C temperature-controlled, 200-ml experimental chambers.29 Three separate cohorts of 20 mice were exposed to 0.90 vol% isoflurane at the same time of day for 2 h on 10 separate occasions over the course of 1 month. This drug concentration was previously determined to be an effective concentration at which > 99% of the population loses the righting reflex under steady-state conditions.17 The righting reflex is a commonly used behavioural endpoint of the anaesthetic state.30 Upon completion of a 2-h anaesthetic challenge, mice were rotated into a supine position, anaesthetic delivery was stopped, and the latency to the return of righting reflex was recorded. Emergence time was defined experimentally as the time when the mouse first spontaneously returned to a prone position. Measurement resolution was to the nearest second. Anaesthetic concentrations were confirmed using a Riken FI-21 refractometer (A.M. Bickford, Wales Center, NY, USA). Normothermia was verified by measuring rectal temperatures. All mice had core temperatures of mean 37.5 (standard deviation [sd] 0.6)°C immediately after anaesthetic emergence.

was to the nearest second. Anaesthetic concentrations were confirmed using a Riken FI-21 refractometer (A.M. Bickford, Wales Center, NY, USA). Normothermia was verified by measuring rectal temperatures. All mice had core temperatures of mean 37.5 (standard deviation [sd] 0.6)°C immediately after anaesthetic emergence. In a previous study, we exposed mice to varying concentrations of isoflurane (0.3, 0.4, 0.6, 0.7 vol%) for 4 h, assessing righting reflex every 3 min (n=20 per concentration).17 The 0.3 and 0.7 vol% isoflurane concentrations were studied with the same 20 mice. We repeated this experiment on four separate occasions.17 Individual response probability was computed as described17,18,31 as the fraction of responsive righting reflex assessments over total number of righting reflex assessments during the final 2 h of the 4-h anaesthetic exposure. Here, we reanalyse these data to determine the stability of anaesthetic sensitivity over experimental days and compute the autocorrelation of individual sensitivity over time.

of responsive righting reflex assessments over total number of righting reflex assessments during the final 2 h of the 4-h anaesthetic exposure. Here, we reanalyse these data to determine the stability of anaesthetic sensitivity over experimental days and compute the autocorrelation of individual sensitivity over time. To identify a probability distribution most likely to explain the observed emergence times, we followed a statistical approach outlined by Clauset and colleagues.32 We considered several commonly used right-skewed distributions as potential models: exponential, lognormal, Weibull, and gamma. The overall goal is to not only determine which of these distributions best fits the data (as there will always be one best fit) but to also determine the likelihood that the best-fit distribution gave rise to the experimental observations. Briefly, the strategy is as follows: (1) estimate the best-fit parameters for a given distribution model of the emergence time data using maximum likelihood methods; (2) quantify the goodness-of-fit using the Kolmogorov–Smirnov (KS) statistic between the empirical and the best-fit distribution; (3) generate random dataset constrained to be the same size as experimental measurements (400 points) by sampling the best-fit distribution and recompute KS (step 2), sampling was repeated 10000 times; and (4) empirically calculate the probability of obtaining the experimentally observed KS by comparing with KS obtained for random datasets.

taset constrained to be the same size as experimental measurements (400 points) by sampling the best-fit distribution and recompute KS (step 2), sampling was repeated 10000 times; and (4) empirically calculate the probability of obtaining the experimentally observed KS by comparing with KS obtained for random datasets. In this case, large P-values indicate experimental observations and a random dataset sampled from the best-fit distribution are statistically indistinguishable. Thus, the goodness-of-fit of experimental data is consistent with finite sample size. To rule in a distribution as a potential model for the experimental observations, we used a threshold of P>0.10.32

erimental observations and a random dataset sampled from the best-fit distribution are statistically indistinguishable. Thus, the goodness-of-fit of experimental data is consistent with finite sample size. To rule in a distribution as a potential model for the experimental observations, we used a threshold of P>0.10.32 A two-compartment PK-PD model12 was constructed to simulate blood and effect-site compartment drug concentrations using the following equations:(1)dc1dt=−(k12+k10)x1+k21c2+I(t)(2)dc2dt=k12x1−k21c2where, c1 and c2 are drug blood and effect-site concentrations, respectively; k12 is the rate constant governing transmission of drug from blood to effect-site; k21 is the drug rate constant from effect-site to blood; k10 is the rate constant of drug elimination; and I(t) is the drug infusion rate, in this case modelled to be a step function lasting 25 time steps chosen such that the first compartment is close to saturation at the end of the exposure. Concentrations and time steps are modelled in arbitrary units. The rate constants used for all PK-PD simulations were: k10=0.8;k12=0.2;k21=0.05.Fig 1Emergence times are highly variable across identical trials in every individual. Return of righting reflex was used as a behavioural measure of anaesthetic emergence for mice exposed to 0.9 vol% isoflurane for 2 h. In total, 40 wild-type male mice were included, all of which were exposed to the same anaesthetic regimen on 10 separate occasions. (a) Overall distribution of emergence times across all mice and trials. (b) Emergence percentile as a function of anaesthetic trial for a single, representative mouse. (c) Emergence times across all mice for all 10 trials. (d) Emergence times across all trials for all 40 mice. (e) Emergence percentile for each trial and mouse combination. (f) The matrix norm of the percentile emergence time matrix (dashed red line) was statistically indistinguishable from shuffled surrogate controls (1000 shuffles, p=0.68, permutation test).Fig 1

trials. (d) Emergence times across all trials for all 40 mice. (e) Emergence percentile for each trial and mouse combination. (f) The matrix norm of the percentile emergence time matrix (dashed red line) was statistically indistinguishable from shuffled surrogate controls (1000 shuffles, p=0.68, permutation test).Fig 1 Emergence times are highly variable across identical trials in every individual. Return of righting reflex was used as a behavioural measure of anaesthetic emergence for mice exposed to 0.9 vol% isoflurane for 2 h. In total, 40 wild-type male mice were included, all of which were exposed to the same anaesthetic regimen on 10 separate occasions. (a) Overall distribution of emergence times across all mice and trials. (b) Emergence percentile as a function of anaesthetic trial for a single, representative mouse. (c) Emergence times across all mice for all 10 trials. (d) Emergence times across all trials for all 40 mice. (e) Emergence percentile for each trial and mouse combination. (f) The matrix norm of the percentile emergence time matrix (dashed red line) was statistically indistinguishable from shuffled surrogate controls (1000 shuffles, p=0.68, permutation test). A PK-derived effect-site concentration was utilised to compute drug effect, E(c2), in the PD component of the model using a standard fixed Hill slope sigmoid concentration–response curve.12(3)E(c2)=1(EC50c2)h

Emergence times are highly variable across identical trials in every individual. Return of righting reflex was used as a behavioural measure of anaesthetic emergence for mice exposed to 0.9 vol% isoflurane for 2 h. In total, 40 wild-type male mice were included, all of which were exposed to the same anaesthetic regimen on 10 separate occasions. (a) Overall distribution of emergence times across all mice and trials. (b) Emergence percentile as a function of anaesthetic trial for a single, representative mouse. (c) Emergence times across all mice for all 10 trials. (d) Emergence times across all trials for all 40 mice. (e) Emergence percentile for each trial and mouse combination. (f) The matrix norm of the percentile emergence time matrix (dashed red line) was statistically indistinguishable from shuffled surrogate controls (1000 shuffles, p=0.68, permutation test). A PK-derived effect-site concentration was utilised to compute drug effect, E(c2), in the PD component of the model using a standard fixed Hill slope sigmoid concentration–response curve.12(3)E(c2)=1(EC50c2)h Here, EC50 is the drug concentration that produces a half-maximal effect, and the Hill slope, h, was set at 15 in accordance with experimental estimates for isoflurane.4,29,33 Given the use of inbred C57BL/6J mice that are nearly genetically identical, we assumed that the PK equations would be identical for all individuals and the variability stems primarily from the PD portion of the model.34 Specifically, E(c2) in equation (3) denotes the probability that an individual will not have righting reflex at a given concentration. Correspondingly, 1−E denotes the probability of having intact righting reflex. This model is the simplest interpretation of the population PD component of the combined PK-PD model, which assumes that the response probability is solely governed by drug concentration in the effect-site compartment. To implement this model at the level of an individual, we simulated a Gaussian random walk (σ=0.1) bounded between 0 and 1 (reflective boundary condition). The points at which this random walk exceeded E, the righting reflex was defined to be present. Thus, emergence time can be defined as the latency from anaesthetic discontinuation until return of the righting reflex. To simulate inter-individual variability, we expressed individual differences in anaesthetic sensitivity as differences in EC50. We simulated 5000 trials (lasting 250 time steps, the first 25 of which are drug administration) at each of the nine different EC50 values from 0.1 to 0.9.

ntil return of the righting reflex. To simulate inter-individual variability, we expressed individual differences in anaesthetic sensitivity as differences in EC50. We simulated 5000 trials (lasting 250 time steps, the first 25 of which are drug administration) at each of the nine different EC50 values from 0.1 to 0.9. As an alternative to the classic PK-PD model based on previous work,12,27 we modelled the PD component as a Markov process containing two behavioural states (‘responsive’ and ‘unresponsive’). Consistent with prior findings at steady-state,17,18 and in a key departure from the classic PD model, each of these states resists transitions. This model relates the current state of a system, X, at time t to that of the previous time step, t−1, as(4)Xt=MXt−1where M is the 2×2 transition probability matrix governing probabilities of state transitions.(5)M=[P(R|R)1−P(R|R)1−P(U|U)P(U|U)]P(R|R) is the conditional probability that a given mouse will stay responsive across two consecutive trials, and P(U|U) is the conditional probability of staying unresponsive across two consecutive trials. Both of these were empirically estimated from the data as described.17 A phenomenon observed in experimental datasets is that the degree of RST is constant for a given anaesthetic across individuals despite large differences in anaesthetic sensitivity.17,18 RST is defined as the average of the main diagonal of M,(6)RST=12∗(P(R|R)+P(U|U))

empirically estimated from the data as described.17 A phenomenon observed in experimental datasets is that the degree of RST is constant for a given anaesthetic across individuals despite large differences in anaesthetic sensitivity.17,18 RST is defined as the average of the main diagonal of M,(6)RST=12∗(P(R|R)+P(U|U)) Conservation of RST constrains M such that knowing a single entry of M will allow for the full description of the transition probability matrix. Solving for P(U|U) in equation (6), we can rewrite M as:(7)M=[P(R|R)1−P(R|R)1−(2∗RST−P(R|R))2∗RST−P(R|R)] Previous studies17,18 described this discrete-time Markov process at steady-state anaesthetic concentrations. However, in the emergence paradigm, by definition, the system will no longer be at steady-state. Stochastic modelling still can be applied to describe the anaesthetic emergence process. Rather than M describing a single, steady-state transition probability matrix, M now represents a family of transition probability matrices as a function of drug brain concentration. Equation (4) can be modified to describe stochastic state switching of emergence from general anaesthesia.(8)Xt=MC(t)Xt−1Here, the current state of system X at time t is related to the state of the system at time t−1, where MC(t) is the transition probability matrix for a given drug brain concentration observed at time t. We assert that at t=0,(9)X0=[01]which denotes that all individuals must start in an unresponsive state.

8)Xt=MC(t)Xt−1Here, the current state of system X at time t is related to the state of the system at time t−1, where MC(t) is the transition probability matrix for a given drug brain concentration observed at time t. We assert that at t=0,(9)X0=[01]which denotes that all individuals must start in an unresponsive state. To generate MC(t), we first needed to model the PK component of emergence from general anaesthesia. We used isoflurane washout kinetics described previously.19 Brain drug concentration was fit to an exponential decay curve to estimate drug washout time,(10)C(t)=C0e−ktwhere, C(t) is drug concentration at time t, C0 is the initial drug brain concentration at the start of the washout period (at t=0), and k is the exponential decay constant. Briefly, mice were exposed to 0.9 vol% isoflurane for 30 min before anaesthetic discontinuation, and subsequent brain harvests were performed at times that were logarithmically spaced for the 1 h after isoflurane cessation. Brain concentrations were determined by high-performance liquid chromatography analysis.19 As in the standard PK-PD model, we assumed that all individuals follow identical washout PK. Best-fit values for C0 and k are shown as:(11)C0=306μgisofluranegbraintissuek=0.39min−1Next, we focused on the PD component of the model. To estimate M across all possible drug concentrations, MC(t), we used experimental estimates of P(R|R) across all drug concentrations. Population-level mean and sd values for P(R|R) were computed across a range of steady-state isoflurane concentrations (0.3, 0.4, 0.6, 0.7, 0.9 vol%).17,18 Two additional experimental conditions were used to complete the simulation. At 0 vol% isoflurane P(R|R) was set to 1. At 1.2 vol% isoflurane P(R|R) was set to 0. These assumptions express the fact that an animal cannot be anaesthetised without anaesthesia and that, at some high enough anaesthetic concentration, all animals are reliably anaesthetised. To obtain MC(t) for all anaesthetic concentrations, we interpolated mean and sd for P(R|R) and RST using a cubic spline. Thus, for each anaesthetic concentration, we can obtain the population mean and the sd for all parameters to completely define MC(t) at each anaesthetic concentration.

als are reliably anaesthetised. To obtain MC(t) for all anaesthetic concentrations, we interpolated mean and sd for P(R|R) and RST using a cubic spline. Thus, for each anaesthetic concentration, we can obtain the population mean and the sd for all parameters to completely define MC(t) at each anaesthetic concentration. To simulate emergence from anaesthesia, at each anaesthetic concentration during the drug washout simulation starting from 0.9 vol% isoflurane, governed by equation (10), we drew a random sample from the distribution of P(R|R) and RST corresponding to the current anaesthetic concentration. This process was repeated 40 times to simulate 40 identical individuals, with 10 runs simulated for each individual. Emergence was defined as the first-time step at which the system arrived in the responsive state.

om sample from the distribution of P(R|R) and RST corresponding to the current anaesthetic concentration. This process was repeated 40 times to simulate 40 identical individuals, with 10 runs simulated for each individual. Emergence was defined as the first-time step at which the system arrived in the responsive state. Data were analysed in MATLAB (2024a, MathWorks, Natick, MA, USA) using the Curve Fitting, Symbolic Math, and Statistics and Machine Learning Toolboxes. Emergence time data were tested for normality using the one-sample KS test. Data were transformed to standardised z-score before normality testing. The χ2 goodness-of-fit test was used to compare individual emergence time distribution with a null hypothesis uniform distribution of population times. Standard quantile bins were used for the expected frequency distribution (bin edges were 0, 25, 50, 75, 100 percentile of population times, expected frequency per bin was 2.5). Correlation analysis of empirical emergence times and anaesthetic sensitivity was computed with Pearson's linear correlation. Autocorrelations of emergence times and anaesthetic sensitivity across experimental days were fit to a single exponential decay curve. Matrix norms from both empirical and model simulated emergence times were compared with matrix norms of randomly shuffled surrogate datasets (1000 shuffles, with replacement). Normalised mutual information was computed using code from Blackwood and colleagues.21 Normalised mutual information was computed for both empirical and model simulated emergence times and compared with results of randomly shuffled surrogate dataset (10 000 shuffles, 30 bins used for joint distribution calculation). P-values were directly computed via permutation test for comparisons with random surrogate shuffles for matrix norm and normalised mutual information comparisons. Other than for empirical emergence time distribution model fitting, P-values <0.05 were considered significant for all comparisons.

Before anaesthetic exposures, mice were habituated to gas-tight, 37°C temperature-controlled, 200-ml experimental chambers.29 Three separate cohorts of 20 mice were exposed to 0.90 vol% isoflurane at the same time of day for 2 h on 10 separate occasions over the course of 1 month. This drug concentration was previously determined to be an effective concentration at which > 99% of the population loses the righting reflex under steady-state conditions.17 The righting reflex is a commonly used behavioural endpoint of the anaesthetic state.30 Upon completion of a 2-h anaesthetic challenge, mice were rotated into a supine position, anaesthetic delivery was stopped, and the latency to the return of righting reflex was recorded. Emergence time was defined experimentally as the time when the mouse first spontaneously returned to a prone position. Measurement resolution was to the nearest second. Anaesthetic concentrations were confirmed using a Riken FI-21 refractometer (A.M. Bickford, Wales Center, NY, USA). Normothermia was verified by measuring rectal temperatures. All mice had core temperatures of mean 37.5 (standard deviation [sd] 0.6)°C immediately after anaesthetic emergence.

In a previous study, we exposed mice to varying concentrations of isoflurane (0.3, 0.4, 0.6, 0.7 vol%) for 4 h, assessing righting reflex every 3 min (n=20 per concentration).17 The 0.3 and 0.7 vol% isoflurane concentrations were studied with the same 20 mice. We repeated this experiment on four separate occasions.17 Individual response probability was computed as described17,18,31 as the fraction of responsive righting reflex assessments over total number of righting reflex assessments during the final 2 h of the 4-h anaesthetic exposure. Here, we reanalyse these data to determine the stability of anaesthetic sensitivity over experimental days and compute the autocorrelation of individual sensitivity over time.

To identify a probability distribution most likely to explain the observed emergence times, we followed a statistical approach outlined by Clauset and colleagues.32 We considered several commonly used right-skewed distributions as potential models: exponential, lognormal, Weibull, and gamma. The overall goal is to not only determine which of these distributions best fits the data (as there will always be one best fit) but to also determine the likelihood that the best-fit distribution gave rise to the experimental observations. Briefly, the strategy is as follows: (1) estimate the best-fit parameters for a given distribution model of the emergence time data using maximum likelihood methods; (2) quantify the goodness-of-fit using the Kolmogorov–Smirnov (KS) statistic between the empirical and the best-fit distribution; (3) generate random dataset constrained to be the same size as experimental measurements (400 points) by sampling the best-fit distribution and recompute KS (step 2), sampling was repeated 10000 times; and (4) empirically calculate the probability of obtaining the experimentally observed KS by comparing with KS obtained for random datasets. In this case, large P-values indicate experimental observations and a random dataset sampled from the best-fit distribution are statistically indistinguishable. Thus, the goodness-of-fit of experimental data is consistent with finite sample size. To rule in a distribution as a potential model for the experimental observations, we used a threshold of P>0.10.32

A two-compartment PK-PD model12 was constructed to simulate blood and effect-site compartment drug concentrations using the following equations:(1)dc1dt=−(k12+k10)x1+k21c2+I(t)(2)dc2dt=k12x1−k21c2where, c1 and c2 are drug blood and effect-site concentrations, respectively; k12 is the rate constant governing transmission of drug from blood to effect-site; k21 is the drug rate constant from effect-site to blood; k10 is the rate constant of drug elimination; and I(t) is the drug infusion rate, in this case modelled to be a step function lasting 25 time steps chosen such that the first compartment is close to saturation at the end of the exposure. Concentrations and time steps are modelled in arbitrary units. The rate constants used for all PK-PD simulations were: k10=0.8;k12=0.2;k21=0.05.Fig 1Emergence times are highly variable across identical trials in every individual. Return of righting reflex was used as a behavioural measure of anaesthetic emergence for mice exposed to 0.9 vol% isoflurane for 2 h. In total, 40 wild-type male mice were included, all of which were exposed to the same anaesthetic regimen on 10 separate occasions. (a) Overall distribution of emergence times across all mice and trials. (b) Emergence percentile as a function of anaesthetic trial for a single, representative mouse. (c) Emergence times across all mice for all 10 trials. (d) Emergence times across all trials for all 40 mice. (e) Emergence percentile for each trial and mouse combination. (f) The matrix norm of the percentile emergence time matrix (dashed red line) was statistically indistinguishable from shuffled surrogate controls (1000 shuffles, p=0.68, permutation test).Fig 1

As an alternative to the classic PK-PD model based on previous work,12,27 we modelled the PD component as a Markov process containing two behavioural states (‘responsive’ and ‘unresponsive’). Consistent with prior findings at steady-state,17,18 and in a key departure from the classic PD model, each of these states resists transitions. This model relates the current state of a system, X, at time t to that of the previous time step, t−1, as(4)Xt=MXt−1where M is the 2×2 transition probability matrix governing probabilities of state transitions.(5)M=[P(R|R)1−P(R|R)1−P(U|U)P(U|U)]P(R|R) is the conditional probability that a given mouse will stay responsive across two consecutive trials, and P(U|U) is the conditional probability of staying unresponsive across two consecutive trials. Both of these were empirically estimated from the data as described.17 A phenomenon observed in experimental datasets is that the degree of RST is constant for a given anaesthetic across individuals despite large differences in anaesthetic sensitivity.17,18 RST is defined as the average of the main diagonal of M,(6)RST=12∗(P(R|R)+P(U|U)) Conservation of RST constrains M such that knowing a single entry of M will allow for the full description of the transition probability matrix. Solving for P(U|U) in equation (6), we can rewrite M as:(7)M=[P(R|R)1−P(R|R)1−(2∗RST−P(R|R))2∗RST−P(R|R)] Previous studies17,18 described this discrete-time Markov process at steady-state anaesthetic concentrations. However, in the emergence paradigm, by definition, the system will no longer be at steady-state.

Conservation of RST constrains M such that knowing a single entry of M will allow for the full description of the transition probability matrix. Solving for P(U|U) in equation (6), we can rewrite M as:(7)M=[P(R|R)1−P(R|R)1−(2∗RST−P(R|R))2∗RST−P(R|R)] Previous studies17,18 described this discrete-time Markov process at steady-state anaesthetic concentrations. However, in the emergence paradigm, by definition, the system will no longer be at steady-state. Stochastic modelling still can be applied to describe the anaesthetic emergence process. Rather than M describing a single, steady-state transition probability matrix, M now represents a family of transition probability matrices as a function of drug brain concentration. Equation (4) can be modified to describe stochastic state switching of emergence from general anaesthesia.(8)Xt=MC(t)Xt−1Here, the current state of system X at time t is related to the state of the system at time t−1, where MC(t) is the transition probability matrix for a given drug brain concentration observed at time t. We assert that at t=0,(9)X0=[01]which denotes that all individuals must start in an unresponsive state.

Data were analysed in MATLAB (2024a, MathWorks, Natick, MA, USA) using the Curve Fitting, Symbolic Math, and Statistics and Machine Learning Toolboxes. Emergence time data were tested for normality using the one-sample KS test. Data were transformed to standardised z-score before normality testing. The χ2 goodness-of-fit test was used to compare individual emergence time distribution with a null hypothesis uniform distribution of population times. Standard quantile bins were used for the expected frequency distribution (bin edges were 0, 25, 50, 75, 100 percentile of population times, expected frequency per bin was 2.5). Correlation analysis of empirical emergence times and anaesthetic sensitivity was computed with Pearson's linear correlation. Autocorrelations of emergence times and anaesthetic sensitivity across experimental days were fit to a single exponential decay curve. Matrix norms from both empirical and model simulated emergence times were compared with matrix norms of randomly shuffled surrogate datasets (1000 shuffles, with replacement). Normalised mutual information was computed using code from Blackwood and colleagues.21 Normalised mutual information was computed for both empirical and model simulated emergence times and compared with results of randomly shuffled surrogate dataset (10 000 shuffles, 30 bins used for joint distribution calculation). P-values were directly computed via permutation test for comparisons with random surrogate shuffles for matrix norm and normalised mutual information comparisons. Other than for empirical emergence time distribution model fitting, P-values <0.05 were considered significant for all comparisons.

We exposed 40 male mice to 0.9 vol% isoflurane for 2 h before measuring emergence time defined as the return of righting reflex, a behavioural correlate of recovery of consciousness in mice. Each mouse was exposed to isoflurane on 10 separate occasions using the same exposure protocol. Overall, the pooled distribution of emergence times was positively skewed (Fig. 1a, skewness >2). The range (9–1387 s) of emergence times from isoflurane spanned more than two orders of magnitude, indicative of a highly variable process. We previously showed17,18 that even in this highly homogeneous population of mice, anaesthetic sensitivity varies significantly. For instance, at a steady-state population EC50 (≈0.6 vol% isoflurane) some individuals were found to be unresponsive on ≈10% of trials, whereas other individuals were unresponsive on ≈90% of trials. Further, individual differences in anaesthetic sensitivity persisted on the time scale of weeks (Fig. 2c and d). If emergence from anaesthesia is primarily determined by drug washout, more sensitive individuals should consistently take a longer time to emerge. In other words, if the variability in emergence times is owing to individual differences in drug sensitivity, one would expect consistent inter-individual differences in emergence times.Fig 2Previous emergence times do not predict future emergence times. (a) There was a weak, positive linear correlation between prior and current emergence times (P=0.043, Pearson's linear correlation). However, the amount of variance explained by correlation was low (R2=0.011). (b) Normalised mutual information for within-subject emergence times was negligible, indicating no statistical dependency on prior emergence times (1000 shuffles, P=0.093, permutation test). (c) There was a strong, positive linear corelation between prior and current estimates of anaesthetic sensitivity across a range of anaesthetic concentrations (0.3 vol% isoflurane [green] to 0.7 vol% isoflurane [blue], R2=0.60, P<0.0001, Pearson's linear correlation).

e times (1000 shuffles, P=0.093, permutation test). (c) There was a strong, positive linear corelation between prior and current estimates of anaesthetic sensitivity across a range of anaesthetic concentrations (0.3 vol% isoflurane [green] to 0.7 vol% isoflurane [blue], R2=0.60, P<0.0001, Pearson's linear correlation). (d) Autocorrelations of anaesthetic sensitivity (red) and emergence times (blue) decay on vastly distinct timescales (sensitivity T1/2=6.90 days (95% CI 5.66–8.84), emergence T1/2=0.66 days (95% CI 0.45–1.24).Fig 2Fig 3Emergence times are highly positively skewed and are best fit by a gamma distribution. (a) Histogram illustrates the distribution of emergence times, indicating a highly positively skewed pattern (skewness =2.08). Potential distributions (lognormal [black], Weibull [red], exponential [blue], and gamma [green]) were fit to the data using maximum likelihood estimation. (b) Fraction of surrogates that fit worse than the experimental data. If <10% of surrogate samples had a KS distance greater than or equal to the dataset, the candidate model was rejected (Plognormal=0.045, PWeibull=0.090, Pexponential=0.043, Pgamma=0.33). KS, Kolmogorov–Smirnov.Fig 3

ing maximum likelihood estimation. (b) Fraction of surrogates that fit worse than the experimental data. If <10% of surrogate samples had a KS distance greater than or equal to the dataset, the candidate model was rejected (Plognormal=0.045, PWeibull=0.090, Pexponential=0.043, Pgamma=0.33). KS, Kolmogorov–Smirnov.Fig 3 Previous emergence times do not predict future emergence times. (a) There was a weak, positive linear correlation between prior and current emergence times (P=0.043, Pearson's linear correlation). However, the amount of variance explained by correlation was low (R2=0.011). (b) Normalised mutual information for within-subject emergence times was negligible, indicating no statistical dependency on prior emergence times (1000 shuffles, P=0.093, permutation test). (c) There was a strong, positive linear corelation between prior and current estimates of anaesthetic sensitivity across a range of anaesthetic concentrations (0.3 vol% isoflurane [green] to 0.7 vol% isoflurane [blue], R2=0.60, P<0.0001, Pearson's linear correlation). (d) Autocorrelations of anaesthetic sensitivity (red) and emergence times (blue) decay on vastly distinct timescales (sensitivity T1/2=6.90 days (95% CI 5.66–8.84), emergence T1/2=0.66 days (95% CI 0.45–1.24).

(0.3 vol% isoflurane [green] to 0.7 vol% isoflurane [blue], R2=0.60, P<0.0001, Pearson's linear correlation). (d) Autocorrelations of anaesthetic sensitivity (red) and emergence times (blue) decay on vastly distinct timescales (sensitivity T1/2=6.90 days (95% CI 5.66–8.84), emergence T1/2=0.66 days (95% CI 0.45–1.24). Emergence times are highly positively skewed and are best fit by a gamma distribution. (a) Histogram illustrates the distribution of emergence times, indicating a highly positively skewed pattern (skewness =2.08). Potential distributions (lognormal [black], Weibull [red], exponential [blue], and gamma [green]) were fit to the data using maximum likelihood estimation. (b) Fraction of surrogates that fit worse than the experimental data. If <10% of surrogate samples had a KS distance greater than or equal to the dataset, the candidate model was rejected (Plognormal=0.045, PWeibull=0.090, Pexponential=0.043, Pgamma=0.33). KS, Kolmogorov–Smirnov.

ta using maximum likelihood estimation. (b) Fraction of surrogates that fit worse than the experimental data. If <10% of surrogate samples had a KS distance greater than or equal to the dataset, the candidate model was rejected (Plognormal=0.045, PWeibull=0.090, Pexponential=0.043, Pgamma=0.33). KS, Kolmogorov–Smirnov. Across all trials, emergence times were highly variable across both individuals (Fig. 1b) and the population (Fig. 1c). To determine whether the overall distribution of emergence times was owing to inter-individual differences, we divided the overall range of emergence times observed over all mice and trials (400 observations) into four quartiles. If one individual consistently emerged sooner than others, we would expect that most of their emergence times would fall into the lower quartiles. Alternatively, if there were no consistent inter-individual differences, we would expect that ∼25% of all observations would fall into each of the four quartiles for each individual. Emergence times in 37 out of 40 mice were statistically indistinguishable (χ2, P>0.05) from uniform distribution among quartiles. This suggests that much of the variability in emergence times is not owing to inter-individual differences.

f all observations would fall into each of the four quartiles for each individual. Emergence times in 37 out of 40 mice were statistically indistinguishable (χ2, P>0.05) from uniform distribution among quartiles. This suggests that much of the variability in emergence times is not owing to inter-individual differences. Repeated exposure to isoflurane did not result in progressive changes in emergence times (Fig. 1d, r=0.0044, P=0.93, Pearson's linear correlation). Variability in emergence times did not differ across individuals or exposures (F(9,39)=0.34, P=0.09, two-sample F-test). This provides further credence to the hypothesis that different anaesthetic exposures are approximately equivalent and that there are no appreciable inter-individual differences in emergence times. To test this hypothesis further, we computed the norm of the matrix of percentile emergence times (Fig. 1e, 40 × 10: animals × trials) and compared it with random shuffles of the same data. If the experimental data that uniquely identifies each individual mouse and each individual trial contains some information about emergence times, then the norm of the matrix ought to be larger than that obtained for shuffle surrogates. Conversely, if experimental observations contain as much information about emergence times as random samples drawn from the overall distribution of emergence times, then the norm of experimental data should be similar to that obtained for random shuffles. The permutation test revealed that the experimentally observed matrix norm (10.38) is statistically indistinguishable from randomly shuffled surrogates (Fig. 1f, P=0.68, permutation test, 1000 shuffles). Altogether, these results suggest that emergence data are statistically indistinguishable from random samples drawn from the overall distribution of emergence times.

ally observed matrix norm (10.38) is statistically indistinguishable from randomly shuffled surrogates (Fig. 1f, P=0.68, permutation test, 1000 shuffles). Altogether, these results suggest that emergence data are statistically indistinguishable from random samples drawn from the overall distribution of emergence times. To determine if there is a sexually dimorphic effect in the observed variability in emergence times,31 we similarly exposed 20 female mice to 10 trials of 0.9 vol% isoflurane for 2 h. The pooled distribution of female emergence times was also highly positively skewed (Supplementary Fig. S1a, skewness >1.5). Compared with male emergence times, the female median emergence time was shorter (Supplementary Fig. S1a, P=0.0013, Mann–Whitney U test). Despite this, female emergence times closely followed the same underlying principal features of high variability at both the population and individual level. Female emergence times varied dramatically over ≈50-fold range, from 21 to 951 s (Supplementary Fig. S1a). Emergence times in 19 out of 20 mice were indistinguishable (Supplementary Fig. S1b, χ2, P>0.05) from uniform distribution among quartiles. Like in the male cohort, repeated exposure to isoflurane did not result in progressive changes in emergence times in females (Supplementary Fig. S1c, r=–0.0052, P=0.94, Pearson's linear correlation). Variability in emergence times did not differ across individuals or exposures either (Supplementary Fig. S1b and c, F(9,19)=0.51, P=0.30, two-sample F-test). Lastly, the observed matrix norm (7.28) of female emergence times (mice × trials) is statistically indistinguishable from randomly shuffled surrogates (Supplementary Fig. S1d and e, P=0.47, permutation test, 1000 shuffles). Together, these results suggest that emergence data are statistically indistinguishable from random samples drawn from the overall distribution of emergence times for both male and female mice.

y indistinguishable from randomly shuffled surrogates (Supplementary Fig. S1d and e, P=0.47, permutation test, 1000 shuffles). Together, these results suggest that emergence data are statistically indistinguishable from random samples drawn from the overall distribution of emergence times for both male and female mice. In standard PK-PD modelling, variability of emergence times is related to differences in anaesthetic sensitivity, drug uptake and elimination, or both. Conversely, in the neuronal dynamics model, emergence times are primarily driven by a quintessentially random process. If we assume, based on previous findings,17 that anaesthetic sensitivity and PK in a given individual change slowly, then the PK-PD model suggests that emergence times observed on two consecutive trials in a given individual ought to be correlated. In contrast, the neural dynamics model predicts that there should not be appreciable correlation. We calculated Pearson correlation coefficients between successive emergence times for each individual across all pairs of consecutive trials (Fig. 2a). Our analysis showed a weak, positive correlation between prior and current emergence times (R2= 0.011, P=0.043, Pearson's linear correlation). This indicates that having information of a prior emergence time would account for only 1% of the variability seen in emergence times. However, because the overall distribution of emergence times is highly right skewed, the P-value for the Pearson correlation might be falsely statistically significant.35

correlation). This indicates that having information of a prior emergence time would account for only 1% of the variability seen in emergence times. However, because the overall distribution of emergence times is highly right skewed, the P-value for the Pearson correlation might be falsely statistically significant.35 Therefore, we implemented a more robust test of inter-trial dependency, normalised mutual information. Unlike correlation analyses, mutual information can detect both linear and nonlinear dependencies and, with appropriate shuffle surrogates, does not strongly depend on the shape of the overall distribution of data. We computed the mutual information between consecutive trials and compared it with mutual information computed on randomly shuffled surrogates (Fig. 2b). This more statistically robust test of inter-trial dependency did not reveal statistical significance (permutation test, P=0.093, compared with 1000 inter-subject shuffled surrogates). Altogether these analyses imply that there is at best a weak relationship between emergence time observed in the same individual on consecutive identical anaesthetic exposures.

ter-trial dependency did not reveal statistical significance (permutation test, P=0.093, compared with 1000 inter-subject shuffled surrogates). Altogether these analyses imply that there is at best a weak relationship between emergence time observed in the same individual on consecutive identical anaesthetic exposures. Variability in emergence times could be related to changes in an individual's anaesthetic sensitivity over time. If fluctuations in anaesthetic sensitivity and emergence times occurred over a similar time scale, it would be challenging to exclude fluctuations in emergence time being dependent on fluctuations in anaesthetic sensitivity. Alternatively, if the rates at which sensitivity and emergence times fluctuate are distinct, anaesthetic sensitivity might not be a primary determinant of anaesthetic emergence. To compare the rates of fluctuations in sensitivity and emergence, we performed serial steady-state righting reflex assessments across multiple concentrations of isoflurane. Mice were subjected to a given anaesthetic concentration on four separate occasions. Our measure of sensitivity was the fraction of trials with an intact righting reflex across each experimental day. In contrast to emergence times, individual sensitivity was strongly positively correlated across sequential anaesthetic exposures across all isoflurane concentrations examined (Fig. 2c, R2=0.60, P<0.0001, Pearson's linear correlation). To determine the stability of correlation across time for both sensitivity and emergence times, we performed an autocorrelation analysis. Autocorrelations of emergence times decayed ∼10 times faster than autocorrelations of sensitivity (Fig. 2d, emergence T1/2=0.66 days [95% CI 0.45–1.24]; sensitivity T1/2=6.90 days [95% CI 5.66–8.84], autocorrelation decay rate). Together, these analyses suggest that the large variability seen in emergence times within individuals cannot be attributed to rapid changes in anaesthetic sensitivity.

tions of sensitivity (Fig. 2d, emergence T1/2=0.66 days [95% CI 0.45–1.24]; sensitivity T1/2=6.90 days [95% CI 5.66–8.84], autocorrelation decay rate). Together, these analyses suggest that the large variability seen in emergence times within individuals cannot be attributed to rapid changes in anaesthetic sensitivity. To analyse how the observed high variability and unpredictability in anaesthetic emergence arise, we pooled emergence times across mice and aimed to find a best-fit distribution model for the emergence dataset (Fig. 3a). We initially found that the data were positively skewed and used a fitting procedure described by Clauset and colleagues.32 We tested several positively skewed distribution models, including the lognormal, Weibull, exponential, and gamma distributions, and found that the gamma distribution model best fit the experimental distribution of emergence times (Fig. 3b, Plognormal=0.045, PWeibull=0.090, Pexponential=0.043, Pgamma=0.33). In this analysis, P>0.10 implies that the gamma distribution fits the data (400 emergence times) and can be expected given the size of the experimental dataset. The only way to identify a better fitting model would be to acquire a significantly larger dataset of emergence times.32Fig 4Standard pharmacokinetic–pharmacodynamic (PK-PD) model of anaesthetic emergence predicts population- but not individual-level features of emergence. (a) Diagram of a two-compartment PK-PD model. I(t) is the drug infusion rate administered into the blood compartment, c1. Drug then can be distributed to the brain, c2, or can be eliminated (c0,notshown). Drug distribution and elimination across compartments are modelled using first-order rate constants, kxy. (b) Blood and brain compartment concentration changes as a function of time using the PK portion of the model after a brief infusion of anaesthetic (25 time steps, arbitrary units [au]). (c) PD component of the model, relating the brain compartment anaesthetic concentration to probability (per unit time) of being anaesthetised. Population-level variability in anaesthetic sensitivity is modelled by a population of shifted dose–response curves (Methods). Blue curve denotes the most sensitive individual, whereas the red curve denotes the most resistant individual. (d) Histogram of simulated emergence time distributions across different anaesthetic sensitivities (coloured histograms). The overall distribution is positively skewed (black line; skewness =2.39).

(Methods). Blue curve denotes the most sensitive individual, whereas the red curve denotes the most resistant individual. (d) Histogram of simulated emergence time distributions across different anaesthetic sensitivities (coloured histograms). The overall distribution is positively skewed (black line; skewness =2.39). (e) In contrast to experimental observations, simulated emergence times are strongly dependent on anaesthetic sensitivity. (f) In contrast to experimental observations, there is a strong positive correlation between simulated prior and current emergence times in each individual (R2=0.88, p<0.0001, Pearson's linear correlation). (g) Normalised mutual information for simulated within-subject (same anaesthetic sensitivity) emergence times was high, indicating strong statistical dependency of current emergence time on prior emergence times (10 000 shuffles, P<0.0001, permutation test). Inset shows normalised mutual information distribution of shuffled surrogate data; red dashed line is simulation normalised mutual information.Fig 1

ivity) emergence times was high, indicating strong statistical dependency of current emergence time on prior emergence times (10 000 shuffles, P<0.0001, permutation test). Inset shows normalised mutual information distribution of shuffled surrogate data; red dashed line is simulation normalised mutual information.Fig 1 Standard pharmacokinetic–pharmacodynamic (PK-PD) model of anaesthetic emergence predicts population- but not individual-level features of emergence. (a) Diagram of a two-compartment PK-PD model. I(t) is the drug infusion rate administered into the blood compartment, c1. Drug then can be distributed to the brain, c2, or can be eliminated (c0,notshown). Drug distribution and elimination across compartments are modelled using first-order rate constants, kxy. (b) Blood and brain compartment concentration changes as a function of time using the PK portion of the model after a brief infusion of anaesthetic (25 time steps, arbitrary units [au]). (c) PD component of the model, relating the brain compartment anaesthetic concentration to probability (per unit time) of being anaesthetised. Population-level variability in anaesthetic sensitivity is modelled by a population of shifted dose–response curves (Methods). Blue curve denotes the most sensitive individual, whereas the red curve denotes the most resistant individual. (d) Histogram of simulated emergence time distributions across different anaesthetic sensitivities (coloured histograms). The overall distribution is positively skewed (black line; skewness =2.39). (e) In contrast to experimental observations, simulated emergence times are strongly dependent on anaesthetic sensitivity. (f) In contrast to experimental observations, there is a strong positive correlation between simulated prior and current emergence times in each individual (R2=0.88, p<0.0001, Pearson's linear correlation). (g) Normalised mutual information for simulated within-subject (same anaesthetic sensitivity) emergence times was high, indicating strong statistical dependency of current emergence time on prior emergence times (10 000 shuffles, P<0.0001, permutation test). Inset shows normalised mutual information distribution of shuffled surrogate data; red dashed line is simulation normalised mutual information.

ensitivity) emergence times was high, indicating strong statistical dependency of current emergence time on prior emergence times (10 000 shuffles, P<0.0001, permutation test). Inset shows normalised mutual information distribution of shuffled surrogate data; red dashed line is simulation normalised mutual information. To quantitatively determine whether a standard PK-PD model could explain the variability in emergence times, we simulated a population of individuals with inter-individual differences in anaesthetic sensitivity. Figure 4a illustrates the PK contribution describing drug distribution across the compartments of the PK-PD model. Drug enters into and is eliminated from the blood compartment, which subsequently distributes drug into the brain effect-site compartment. Blood and brain drug concentration change as a function of time and infusion rate, governed by a system of first-order differential equations (Fig. 4b). The PD portion of the model relates drug effect, the probability of being anaesthetised, to drug brain concentration using a standard Hill equation curve (Fig. 4c). Individual variability in anaesthetic sensitivity was accounted for in model simulations through the shift of the EC50 of the dose–response curve.17,34 This model was then used to simulate the emergence process. The standard PK-PD model was able to recapitulate the overall right-skewed population-level distribution of emergence times (Fig. 4d, skewness >2). However, it failed to capture several salient features of emergence times at a level of an individual. In contrast to experimental observations, the standard PK-PD model predicts that inter-trial variability within an individual should be smaller than inter-individual differences (F(8,44991)=0.127070, P<0.0001, Brown–Forsythe test). Further, the PK-PD model predicts that the average (across trial) emergence time is strongly dependent on the individual (Fig. 4e, F(8,44991)=82002, P<0.0001, one-way anova) and that prior and future emergence times are strongly positively correlated (Fig. 4f, R2=0.88, p<0.0001, Pearson's linear correlation). Lastly, standard PK-PD modelling of emergence predicts significant within-subject mutual information across emergence times (Fig. 4g, P<0.0001, permutation test). Cumulatively, these results demonstrate standard PK-PD modelling reproduces population-level features of emergence but yields incorrect predictions at the level of individuals.

-PD modelling of emergence predicts significant within-subject mutual information across emergence times (Fig. 4g, P<0.0001, permutation test). Cumulatively, these results demonstrate standard PK-PD modelling reproduces population-level features of emergence but yields incorrect predictions at the level of individuals. Neuronal dynamics represented as diffusion on a continuous energy landscape with two potential wells signifying the ‘wake/responsive’ and ‘anaesthetised/unresponsive’ states has been proposed,12,27 and demonstrated17,18 as an effective model for explaining neural inertia and RST. However, this model has only been experimentally verified at steady-state drug concentrations. We aimed to extend this model to describe the nonequilibrium behaviour during emergence from anaesthesia.

sive’ states has been proposed,12,27 and demonstrated17,18 as an effective model for explaining neural inertia and RST. However, this model has only been experimentally verified at steady-state drug concentrations. We aimed to extend this model to describe the nonequilibrium behaviour during emergence from anaesthesia. We constructed a model describing a two-state Markov process that incorporates nonequilibrium PK and conservation of RST into the transition probability matrix governing the likelihood of transitioning between the unresponsive and responsive states (Fig. 5a). Brain drug concentrations at any given time were estimated using isoflurane washout kinetics fitted to an exponential decay function.19 For each timepoint in which the Markov process is evaluated and the transition probability matrix used to estimate the system state was selected based on the corresponding drug concentration at that time (Fig. 5b). Consistent with experimental observations17,18 and the standard PK-PD simulation, this model also included individual differences in anaesthetic sensitivity (see Methods). The only salient difference between the PK-PD and the neuronal dynamics model is that the PK-PD model assumes that drug concentration directly determines the probability of response, whereas the neuronal dynamics model includes RST.Fig 5Neuronal dynamics model of anaesthetic emergence predicts both population- and individual-level features of emergence. (a) Schematic illustrating a bistable neuronal dynamics model. A network consists of two populations of neurones that are self-excitatory and mutually inhibitory. We define these populations as either being active during the unresponsive state (blue) or active during the responsive state (green). Activity of this bistable model can be parametrised as an energy function (top), where well depths are directly related to the probability of observing an animal in either the responsive (R) or unresponsive (U) state. Anaesthetic concentration is shown by shading from blue (high concentration) to green (low concentration). The shape of this energy function depends on the drug concentration. In the limit of high concentration only the unresponsive state is stable. In the limit of low concentration only the responsive state is stable. At the intermediate concentration range, both U and R have comparable stability. Dynamics on this energy landscape were approximated by a transition probability matrix (Methods).

the limit of high concentration only the unresponsive state is stable. In the limit of low concentration only the responsive state is stable. At the intermediate concentration range, both U and R have comparable stability. Dynamics on this energy landscape were approximated by a transition probability matrix (Methods). Changes in the shape of the energy landscape were driven by exponential decay of isoflurane estimated from data. (b) Transition probabilities of staying in the responsive state (P(R|R) (green) or unresponsive state P(U|U) (blue) as a function of time across 40 simulated individuals. (c) Overall simulated emergence time distributions recapitulate positive skew (skewness =1.00). (d) Emergence times for 40 simulations including 10 trials for each simulated mouse. Critically, the family of transition probability matrices (one for each concentration) were identical on all trials for each simulated mouse. Thus, the observed inter-trial variability in emergence times is attributed solely to the stochastic nature of the model. (e) The matrix norm of the emergence time matrix is statistically indistinguishable from random shuffled surrogate controls as seen in the empirical data (1000 shuffles, p=0.56, permutation test).Fig 5

s, the observed inter-trial variability in emergence times is attributed solely to the stochastic nature of the model. (e) The matrix norm of the emergence time matrix is statistically indistinguishable from random shuffled surrogate controls as seen in the empirical data (1000 shuffles, p=0.56, permutation test).Fig 5 Neuronal dynamics model of anaesthetic emergence predicts both population- and individual-level features of emergence. (a) Schematic illustrating a bistable neuronal dynamics model. A network consists of two populations of neurones that are self-excitatory and mutually inhibitory. We define these populations as either being active during the unresponsive state (blue) or active during the responsive state (green). Activity of this bistable model can be parametrised as an energy function (top), where well depths are directly related to the probability of observing an animal in either the responsive (R) or unresponsive (U) state. Anaesthetic concentration is shown by shading from blue (high concentration) to green (low concentration). The shape of this energy function depends on the drug concentration. In the limit of high concentration only the unresponsive state is stable. In the limit of low concentration only the responsive state is stable. At the intermediate concentration range, both U and R have comparable stability. Dynamics on this energy landscape were approximated by a transition probability matrix (Methods). Changes in the shape of the energy landscape were driven by exponential decay of isoflurane estimated from data. (b) Transition probabilities of staying in the responsive state (P(R|R) (green) or unresponsive state P(U|U) (blue) as a function of time across 40 simulated individuals. (c) Overall simulated emergence time distributions recapitulate positive skew (skewness =1.00). (d) Emergence times for 40 simulations including 10 trials for each simulated mouse. Critically, the family of transition probability matrices (one for each concentration) were identical on all trials for each simulated mouse. Thus, the observed inter-trial variability in emergence times is attributed solely to the stochastic nature of the model. (e) The matrix norm of the emergence time matrix is statistically indistinguishable from random shuffled surrogate controls as seen in the empirical data (1000 shuffles, p=0.56, permutation test).

. Thus, the observed inter-trial variability in emergence times is attributed solely to the stochastic nature of the model. (e) The matrix norm of the emergence time matrix is statistically indistinguishable from random shuffled surrogate controls as seen in the empirical data (1000 shuffles, p=0.56, permutation test). Using the neuronal dynamics model, we simulated emergence from anaesthesia in 40 individuals repeatedly exposed to isoflurane (10 times). Consistent with experimental observations, the population-level distribution of emergence times was positively skewed (Fig. 5c, skewness =1.0). However, unlike the standard PK-PD model, individual-level features of emergence are fully recapitulated using neuronal dynamics modelling. Consistent with experimental observations, we did not observe inter-individual differences in emergence times (Fig. 5d). The distribution of emergence times was statistically indistinguishable from a randomly generated process (Fig. 5e, P=0.56, permutation test). These results collectively demonstrate that neuronal dynamics modelling accurately replicate both population- and individual-level features of anaesthetic emergence.

To quantitatively determine whether a standard PK-PD model could explain the variability in emergence times, we simulated a population of individuals with inter-individual differences in anaesthetic sensitivity. Figure 4a illustrates the PK contribution describing drug distribution across the compartments of the PK-PD model. Drug enters into and is eliminated from the blood compartment, which subsequently distributes drug into the brain effect-site compartment. Blood and brain drug concentration change as a function of time and infusion rate, governed by a system of first-order differential equations (Fig. 4b). The PD portion of the model relates drug effect, the probability of being anaesthetised, to drug brain concentration using a standard Hill equation curve (Fig. 4c). Individual variability in anaesthetic sensitivity was accounted for in model simulations through the shift of the EC50 of the dose–response curve.17,34 This model was then used to simulate the emergence process. The standard PK-PD model was able to recapitulate the overall right-skewed population-level distribution of emergence times (Fig. 4d, skewness >2). However, it failed to capture several salient features of emergence times at a level of an individual. In contrast to experimental observations, the standard PK-PD model predicts that inter-trial variability within an individual should be smaller than inter-individual differences (F(8,44991)=0.127070, P<0.0001, Brown–Forsythe test). Further, the PK-PD model predicts that the average (across trial) emergence time is strongly dependent on the individual (Fig. 4e, F(8,44991)=82002, P<0.0001, one-way anova) and that prior and future emergence times are strongly positively correlated (Fig. 4f, R2=0.88, p<0.0001, Pearson's linear correlation). Lastly, standard PK-PD modelling of emergence predicts significant within-subject mutual information across emergence times (Fig. 4g, P<0.0001, permutation test). Cumulatively, these results demonstrate standard PK-PD modelling reproduces population-level features of emergence but yields incorrect predictions at the level of individuals.

Neuronal dynamics represented as diffusion on a continuous energy landscape with two potential wells signifying the ‘wake/responsive’ and ‘anaesthetised/unresponsive’ states has been proposed,12,27 and demonstrated17,18 as an effective model for explaining neural inertia and RST. However, this model has only been experimentally verified at steady-state drug concentrations. We aimed to extend this model to describe the nonequilibrium behaviour during emergence from anaesthesia.

Emergence times in individuals exposed repeatedly to an identical anaesthetic regimen varied over two orders of magnitude. Although our previous work17,18,31 identified large individual differences in anaesthetic sensitivity, which are classically predicted to be associated with individual differences in emergence times, we failed to observe such individual differences in emergence times. Rather, to a good approximation, emergence times were found to be random samples from the gamma distribution. These observations are a stark departure from the standard PK-PD formulation of anaesthetic effects. Although the PK-PD model recapitulates the behaviour of a population, it yields incorrect predictions at an individual level. In contrast, a neuronal dynamics model that allows for the possibility of neural inertia readily recapitulates both the population- and individual-level emergence time data. Cumulatively, our results imply that emergence from anaesthesia is driven by a quintessentially stochastic process.

ctions at an individual level. In contrast, a neuronal dynamics model that allows for the possibility of neural inertia readily recapitulates both the population- and individual-level emergence time data. Cumulatively, our results imply that emergence from anaesthesia is driven by a quintessentially stochastic process. Standard PK-PD models suggest that deeper anaesthesia should lead to longer recovery times. However, data from humans show that individuals anaesthetised so deeply that they exhibit EEG burst suppression emerge as quickly as those less deeply anaesthetised.36 This discrepancy might be explained by RST, a hidden variable not accounted for in population-level pharmacology. RST is the key property of the neuronal dynamics model that distinguishes it from the standard PK-PD model. Recent experiments at anaesthetic steady-state showed that fluctuations in responsiveness exhibit inertia: an individual is more likely to remain in its previously determined state than undergo a state change.17,18 We hypothesised that this inertia should increase the variability in the time of return to consciousness across identical exposures to anaesthesia.27 Our results showing that emergence times span two orders of magnitude in a single individual are consistent with this hypothesis. The neuronal dynamics model predicts that individual sensitivity and emergence times should be largely independent, and we demonstrate that large fluctuations in emergence times cannot be accounted for by similar changes in anaesthetic sensitivity. Prior examinations of burst suppression as a measure of anaesthetic sensitivity corroborate the lack of intraoperative measures of anaesthetic depth with time to emergence from anaesthesia.31,36

demonstrate that large fluctuations in emergence times cannot be accounted for by similar changes in anaesthetic sensitivity. Prior examinations of burst suppression as a measure of anaesthetic sensitivity corroborate the lack of intraoperative measures of anaesthetic depth with time to emergence from anaesthesia.31,36 A key observation is that emergence times appear to be randomly distributed across individuals and trials, aligning with the concept of sampling from a maximum entropy probability distribution. Notably, the overall population distribution for emergence times is best fit by a gamma distribution, a maximum entropy model that requires minimal prior information for its description.37 Clinical studies examining distributions of emergence times report comparable heavy-tailed distribution fits.38,39 This observation is significant as gamma distributions often arise in processes governing waiting or leaving times and can result from the sum of exponential distributions.40 Our nonequilibrium neuronal dynamics model can be thought of as the sum of individual leaving rates from a series of 2×2 transition probability matrices. In this model, leaving rates for each row of these matrices follow an exponential distribution, leading to a prediction of emergence times that follow a gamma distribution. This result suggests that emergence from anaesthesia is ergodic, meaning that a distribution of many emergence times from a single individual would be equivalent to the combined distribution of a single emergence time from many individuals. Ergodicity is a fundamental concept of stochastic processes,40 and we observe this phenomenon in our collected emergence times, where the emergence time distributions of single individuals mirror those of single trials across many individuals. A key prediction for future studies will be to test whether the ratio of measures of spread, such as variance, between individual and trial emergence time distributions is unitary.

lected emergence times, where the emergence time distributions of single individuals mirror those of single trials across many individuals. A key prediction for future studies will be to test whether the ratio of measures of spread, such as variance, between individual and trial emergence time distributions is unitary. Determining whether a process is truly ergodic requires sufficiently large number of measurements. Although we did not perform a priori power analyses to determine the exact number of observations needed, several factors worked in our favour for observing this effect, including a fast convergence rate (i.e. observations of a system that quickly transitions from unresponsive to responsive) and a low dimensional state space with relatively simple dynamics. Additionally, we used numerous statistically rigorous tests that took finite sample size into account. Despite our best efforts, it remains a possibility that increasing the experimental sample size could result in changes to our best-fit emergence time gamma distribution model. Future studies examining emergence times as a critical measure should be mindful that low sample sizes could lead to incorrect results and conclusions. Ensuring a sufficiently large sample size is essential to capture the dynamics of anaesthetic emergence accurately and avoid potential misinterpretations.

bution model. Future studies examining emergence times as a critical measure should be mindful that low sample sizes could lead to incorrect results and conclusions. Ensuring a sufficiently large sample size is essential to capture the dynamics of anaesthetic emergence accurately and avoid potential misinterpretations. Our neuronal dynamics model describing the nonequilibrium PK state during emergence from general anaesthesia provides a natural extension of previously proposed models of stochastic state switching.12,17,18,27 How well this model can describe the induction process remains unclear. Compared with washout, the wash-in kinetics of isoflurane are faster,19 but this could be accounted for in the model. Something that is less easily accounted for is the PD component. We measure emergence as the time it takes the mouse to turn from supine to prone, or the spontaneous return of righting reflex. To measure induction as the time to the first loss of righting reflex, the study design would require continuous righting reflex assessments. Continuous righting reflex assessments, however, might introduce a confound by artificially increasing the conditional probability of being responsive and biasing the system towards the responsive state. This limitation of our current approach could be mitigated by use of other behavioural and neurophysiological measures.30,41 Future studies are needed to explore stochastic state switching models on the induction process of anaesthetic-induced unconsciousness.

onsive and biasing the system towards the responsive state. This limitation of our current approach could be mitigated by use of other behavioural and neurophysiological measures.30,41 Future studies are needed to explore stochastic state switching models on the induction process of anaesthetic-induced unconsciousness. Emergence from anaesthesia is a complex process influenced by numerous external factors. Variables such as the type of anaesthetic agent,42 age,43 sex,31 stress,44 sleep deprivation,45 pre-existing health conditions,46 genetic variability,5 metabolic interactions,47 and hypothermia48 can all impact emergence time. Our study was specifically designed to minimise independent variables and examine emergence under carefully controlled conditions. Despite this, anaesthetic emergence times still varied over two orders of magnitude, even within the same individual. Future research should target the underlying mechanisms of variability in emergence, even in controlled settings. Insights from such studies could be applied in clinical practice to improve the reliability and predictability of the emergence process.

varied over two orders of magnitude, even within the same individual. Future research should target the underlying mechanisms of variability in emergence, even in controlled settings. Insights from such studies could be applied in clinical practice to improve the reliability and predictability of the emergence process. Our work offers a novel framework that illuminates a root cause of delayed emergence. Specifically, we showed that the causes of delayed emergence are not related to individual differences but are related to the fluctuating internal state of the animal. We described this internal state fluctuation as a random process because, to a good approximation, we failed to identify dependence between individual emergence times. This does not, however, mean that the underlying causes of such state fluctuations are inherently random. Instead, it is possible that state fluctuations under anaesthesia are related to endogenous fluctuations in arousal. Consistent with this conjecture, modulation of specific nodes in this vast arousal network accelerate49, 50, 51, 52, 53, 54, 55 or prolong5,56, 57, 58, 59 emergence from anaesthesia. Because the arousal circuit in the brain is a vast interconnected network composed of many neurones, at the macroscopic level of behavioural observations its effect is indistinguishable from that of an uncorrelated stochastic process. Detailed investigations of the neuronal mechanisms that give rise to such apparently stochastic fluctuations in emergence time might lead to therapeutic interventions aimed at making emergence from general anaesthesia more reliable.

10.13039/100030719National Institute of Health grants (T32 HL007953 to MES, R01 GM144377 to MBK/AP, R01 GM151556 to AP/MBK, and T32 GM112596-09 to AZW); Department of Anesthesiology and Critical Care at the University of Pennsylvania.