Seven domestic pigs (weight 31.1 SD 1.6 kg, 6 male, 1 female) were studied under general anaesthesia and muscle relaxation at the Hedenstierna Laboratory, University of Uppsala, Sweden. Ethical approval was granted by the Uppsala Regional Animal Research Ethics Committee (ref. C98/16). The relevant sections of the ARRIVE guidelines [11] were adhered to.

Experimental protocol

Following premedication with intramuscular xylazine 2 mg/kg, ketamine 20 mg/kg and midazolam 0.5 mg/kg, an ear vein was cannulated, and general anaesthesia commenced and maintained with an infusion of ketamine 32 mg/kg/h, fentanyl 4 mcg/kg/h and midazolam 0.16 mg/kg/h. Muscle relaxation was maintained with an intravenous infusion of rocuronium 30–100 mg/h titrated against spontaneous respiratory effort. Adequacy of anaesthesia was confirmed by absence of reaction to painful stimulation between the front hooves, and absence of any signs of sympathetic stimulation after administration of rocuronium. Aliquots of fentanyl 100–300 mcg were administered as required.

Following induction of general anaesthesia, an endotracheal tube was placed via tracheostomy and mechanical ventilation commenced using a Servo-I ventilator (Maquet, Rastatt, Germany). Immediately following tracheostomy, to maintain normoxia and normocapnia, the ventilation settings were fixed to volume control mode, tidal volume 8 mL/kg, respiratory rate 25/min, PEEP 5 cmH2O and FiO2 sufficient to maintain peripheral oxygen saturations above 92%. Settings were subsequently adjusted for the recording of PaO2 oscillations (see below). Ringer’s lactate solution was infused at 20 mL/kg/h for the first hour then reduced to 10 mL/kg/h for the remainder of the experiment.

These experiments were performed on a subset of animals included in a previously reported study [12] in which the indwelling fast responding PaO2 sensor was available. In this previous study, animals were randomised to either the lung-injury or uninjured group using simple randomisation. In the lung-injury group, the saline-lavage surfactant depletion technique of Lachmann [13] was used aiming for a PaO2/FiO2 ratio of 20–26.7 kPa (150–200 mmHg). 30 mL/kg of warmed saline was administered intratracheally until the PaO2 fell to less than 6.6 kPa (50 mmHg) or mean pulmonary artery pressure rose above 40 mmHg, at which point the lavage fluid was drained from the animal’s lungs and mechanical ventilation was restarted. The lavage process was repeated if PaO2 subsequently rose above 26.7 kPa (200 mmHg). PaO2 was recorded prior to each lavage, immediately after restarting ventilation and subsequently at 30 s intervals. FiO2 was fixed at 1.0 for the entire lavage procedure. The currently reported experiments were performed after induction of lung injury (if applicable), but prior to transfer to the CT scanner, where the previously reported study [12] was performed. Aside from presence/absence of experimental lung injury, there was no section of the protocol for either set of experiments that would adversely affect the results of the other study.

Each animal was ventilated at 12 different ventilatory conditions, representing all combinations of PEEP 5, 8, 10 and 12 cmH2O and tidal volume 7, 10 and 15 mL/kg. For each animal the order of ventilatory conditions was randomised using simple randomisation on pre-printed case record forms. Given the exploratory nature of this study the researchers were not blinded to group (lung-injury vs uninjured) allocation or the order of ventilatory conditions. Each animal was ventilated at each condition with respiratory rate 10/min, inspiratory:expiratory ratio 1:2 in pressure control mode with FiO2 sufficient to maintain PaO2 between 13.3 and 66.5 kPa (100 to 500 mmHg). This lower limit was chosen to ensure haemoglobin was fully saturated for the duration of the experiments so any intra-tidal variation in arterial blood oxygen content was directly proportional to PaO2. If this constraint is not observed, then the PaO2 waveform is non-linearly damped, because any intra-tidal variation in arterial oxygen content would to some extent be buffered by O2 release/absorption by haemoglobin [14]. The upper PaO2 limit of 66.5 kPa (500 mmHg) represented the upper limit of the PaO2 sensor used.

This was a pilot study to investigate the technique used to identify different PaO2 oscillation morphologies. Therefore, sample size was chosen to be sufficient to generate multiple morphologies of PaO2 oscillation rather than to demonstrate statistically significant differences between any identified PaO2 morphology. Based upon the results of previous experiments [8, 9] using fast-acting intra-arterial sensors to detect PaO2 oscillations, between 5 and 8 pigs would be sufficient to demonstrate a wide variety of PaO2 oscillation morphologies.

Measurements

The left common carotid artery was cannulated using a standard 20G 80 mm arterial cannula (LeaderCath, Vygon, Ecouen, France). Through this cannula, a fibreoptic PaO2 sensor with a response time of < 100 ms was introduced and interfaced to an OxyLite Pro (Oxford Optronix, Oxford, UK). These sensors have been described previously [15, 16]. The analogue output from these sensors was continuously recorded (PowerLab and LabChart, ADInstruments, Dunedin, New Zealand) at a sampling frequency of 10 Hz.

Systemic arterial blood pressure was measured using a catheter placed in a branch of the right external carotid artery. A central venous catheter was sited via the right internal jugular vein to measure central venous pressure. A pulmonary artery catheter (Criticath SP5107U, Merit Medical, South Jordan, UT, USA) was placed via the right internal jugular vein using standard flotation methods. All blood pressure waveforms, as well as ECG and tail pulse oximetry, were continuously recorded (IntelliVue M8004A, Philips Healthcare, Best, Netherlands). Cardiac output was measured using a standard thermodilution bolus technique.

Respiratory waveforms including pressure and flow (with tidal volume calculated from the latter) were recorded using a Capnomac Ultima (Datex-Ohmeda, Madison, WI). Mechanical power was calculated using the simplified method for pressure control ventilation [17].

Alignment of PaO2 oscillations with the phases of ventilation

For each ventilatory condition, the oscillations were aligned with the phases of ventilation (inspiration and expiration) by studying the PaO2 increase following the resumption of ventilation after a 20 s breath hold at end-expiration. The period between the start of inspiration (as measured by the airway pressure trace) and the first increase in PaO2, after the steady decline during the breath hold, was measured (Additional file 1: Fig. S1). This period (offset) was then retrospectively applied to the time axis of the PaO2 signal during mechanical ventilation immediately before the breath hold.

Criteria for including PaO2 oscillations in analysis

PaO2 signals from certain ventilatory conditions were excluded from further analyses based upon pre-defined criteria. Signals with trough values less than 13.3 kPa (100 mmHg) were excluded to avoid non-linear dampening of the signal by O2 buffering by haemoglobin as described above. Others were excluded due to the signal-to-noise ratio of the PaO2 oscillations being less than 20 dB (empirically chosen based on preliminary results as being the minimum required to adequately detect the respiratory rate from the signal in 90% of cases). Finally, those oscillations where it was impossible to align the PaO2 signal with the phases of ventilation were excluded.

PaO2 oscillation analysis

The PaO2 oscillations recorded for each ventilatory condition were assigned into one of several arbitrary clusters based upon the morphology of an aggregate single tidal breath. PaO2 oscillations from ten complete breaths were extracted for each ventilatory condition. Any linear trend in the ten-breath signal was removed. The ten breaths were combined into a single breath using a median of the PaO2 of each individual timepoint within the breaths to remove changes in PaO2 occurring at other frequencies, for example, the heart rate. This single breath was then normalized by converting the raw amplitude values (in mmHg) into multiples of the standard deviation of that breath. This conversion allowed comparison of the shape of the oscillations between different ventilatory conditions, without these being impacted by the amplitude of the oscillation or absolute PaO2 value.

Functional principal component analysis

Functional principal component analysis (FPCA) [18] was then used to identify common signals in these normalised breaths. FPCA is a method of dimensionality reduction for complex time-series data, where the mean input PaO2 waveform is calculated and subtracted from each input waveform. A series of eigenfunctions is then created, such that all input waveforms can be reliably described by a weighted sum of these eigenfunctions. The mean waveform and eigenfunctions are constant between all input PaO2 waveforms, but the weights, termed principal component estimates, vary. Thus, a complex input waveform with 600 data points can be described by only a handful of principal components. The FPCA algorithm attempts to choose the optimum eigenfunctions, such that the maximum variance in the input waveforms can be explained with the fewest principal components. The fdapace [19] package in the R statistical software version 4.1.2 (R Foundation for Statistical Computing, Vienna, Austria) was used to calculate functional principal components.

k-means clustering

The principal component estimates for each ventilatory condition were used to assign each ventilatory condition to a cluster using k-means clustering [20]. This attempts to partition each set of inputs (here the principal components) into unique clusters, thus each PaO2 input waveform is assigned to a single cluster. In our data, each point is multi-dimensional comprising a certain number of principal components. k-means clustering generates a certain number of centre points, such that, for each cluster, the sum of the squared distances between each set of input principal component estimates and its cluster centre is minimized. We used the iterative k-means clustering algorithm of Hartigan and Wong [21] to assign individual PaO2 waveforms to clusters.

Choice of both the number of input dimensions (number of principal components) as well as the number of clusters is problematic in k-means clustering as both can significantly affect the quality of the output. We initially set both the number of principal components and number of clusters as the number of principal components required to explain ≥ 95% of the variance in the input waveforms (semantically equivalent to a R2 value of ≥ 0.95 in other statistical models) and then optimised these based upon the Akaike Information Criterion (AIC) [22]. AIC was calculated based upon k-means clustering outputs [23]. Given that k-means clustering is stochastic, Monte Carlo simulation was performed over 10,000 iterations to generate AIC estimates for all potential cluster numbers ranging from zero to the maximum number of principal components identified.

Statistical analyses

Differences in mean PaO2, PaO2 oscillation peak-to-trough height, presence of lung injury and cardiorespiratory variables between the identified PaO2 clusters were analysed using one-way analysis of variance or chi-squared testing as appropriate. Post-hoc testing was not performed between the individual clusters as this was a pilot study and inadequately powered to identify any such differences. Following visual allocation of cluster morphologies into a potential progression, we undertook linear regression analysis of the effects of this potential progression upon principal component estimates. In addition, an exploratory analysis of the effect of any unmeasured intra-animal effect upon cluster membership was performed. A mixed-effects linear model was developed with potential cluster progression as a continuous dependent variable, animal as a random-effects independent variable and other measured baseline cardiorespiratory parameters as fixed-effects dependent variables. The selection of fixed-effects variables was chosen, such that AIC was optimised. For these data, AIC was optimised when PaO2/FiO2 ratio and pulse pressure variation were included as fixed effects variables. Analyses were performed using the R statistical software version 4.1.2 (R Foundation for Statistical Computing, Vienna, Austria).