Digital Twins of Acute Hypoxemic Respiratory Failure Patients Suggest a Mechanistic Basis for Success and Failure of Noninvasive Ventilation

OBJECTIVES: To clarify the mechanistic basis for the success or failure of noninvasive ventilation (NIV) in acute hypoxemic respiratory failure (AHRF). DESIGN: We created digital twins based on mechanistic computational models of individual patients with AHRF. SETTING: Interdisciplinary Collaboration in Systems Medicine Research Network. SUBJECTS: We used individual patient data from 30 moderate-to-severe AHRF patients who had failed high-flow nasal cannula (HFNC) therapy and subsequently underwent a trial of NIV. INTERVENTIONS: Using the digital twins, we evaluated lung mechanics, quantified the separate contributions of external support and patient respiratory effort to lung injury indices, and investigated their relative impact on NIV success or failure. MEASUREMENTS AND MAIN RESULTS: In digital twins of patients who successfully completed/failed NIV, after 2 hours of the trial the mean (sd) of the change in total lung stress was –10.9 (6.2)/–0.35 (3.38) cm H2O, mechanical power –13.4 (12.2)/–1.0 (5.4) J/min, and total lung strain 0.02 (0.24)/0.16 (0.30). In the digital twins, positive end-expiratory pressure (PEEP) produced by HFNC was similar to that set during NIV. In digital twins of patients who failed NIV vs. those who succeeded, intrinsic PEEP was 3.5 (0.6) vs. 2.3 (0.8) cm H2O, inspiratory pressure support was 8.3 (5.9) vs. 22.3 (7.2) cm H2O, and tidal volume was 10.9 (1.2) vs. 9.4 (1.8) mL/kg. In digital twins, successful NIV increased respiratory system compliance +25.0 (16.4) mL/cm H2O, lowered inspiratory muscle pressure –9.7 (9.6) cm H2O, and reduced the contribution of patient spontaneous breathing to total driving pressure by 57.0%. CONCLUSIONS: In digital twins of AHRF patients, successful NIV improved lung mechanics, lowering respiratory effort and indices associated with lung injury. NIV failed in patients for whom only low levels of positive inspiratory pressure support could be applied without risking patient self-inflicted lung injury due to excessive tidal volumes.


Table of Contents
The model employed in this paper has been developed over the past several years and has been applied and validated on a number of different studies [E1]-[E10].The model is organised as a system of several components (see Fig S1 .1), each component representing different aspects of pulmonary dynamics and blood gas transport, e.g. the transport of air in the mouth, the tidal flow in the airways, the gas exchange in the alveolar compartments and their corresponding capillary compartment, the flow of blood in the arteries, the veins, the cardiovascular system, and the gas exchange process in the peripheral tissue compartments.Each component is described as several mass conserving functions and solved as algebraic equations, obtained or approximated from the published literature, experimental data and clinical observations.These equations are solved in series in an iterative manner, so that solving one equation at current time instant (t k ) determines the values of the independent variables in the next equation.At the end of the iteration, the results of the solution of the final equations determine the independent variables of the first equation for the next iteration.
The iterative process continues for a predetermined time, T, representing the total simulation time, with each iteration representing a 'time slice' t of real physiological time (set to 10 ms).At the first iteration(t k , k = 0), an initial set of independent variables are chosen based on values selected by the user.The user can alter these initial variables to investigate the response of the model or to simulate different pathophysiological conditions.Subsequent iterations (t k = t k−1 + ) update the model parameters based on the equations below.
The pulmonary model consists of the mechanical ventilation equipment, anatomical and alveolar dead space, anatomical and alveolar shunts, ventilated alveolar compartments and corresponding perfused capillary compartments.The pressure differential created by the mechanical ventilator or inspiratory muscles (i.e. when modelling spontaneous breathing) drives the flow of gas through the system.The series dead space (SD) is located between the mouth and the alveolar compartments and consists of the trachea, bronchi and the bronchioles where no gas exchange occurs.Inhaled gases pass through the SD during inspiration and alveolar gases pass through the SD during expiration.In the model, an SD of volume 60ml is split into 50 stacked layers of equal volumes (  = 50).No mixing between the compartments of the SD is assumed.
Any residual alveolar air in the SD at the end of expiration is re-inhaled as inspiration is initiated.This residual air is composed of gases exhaled from both perfused alveolar compartments (normal perfusion) and the parallel dead space (PD) (alveolar compartments with limited perfusion).Therefore, the size of dead space (SD and PD) can have a significant effect on the gas composition of the alveolar compartments.
The inhaled air is initially assumed to consist of five gases: oxygen (O2), nitrogen (N2), carbon dioxide (CO2), water vapor (H2O) and a 5th gas (α) used to model additives such as helium or other anaesthetic gases.During an iteration of the model, the flow (f) of air to or from an alveolar compartment i at time t k is determined by the following equation: where   (t k ) is the pressure at the mouth at (t k ),   (t k ) is the pressure in the alveolar compartment  at (t k ), Ru is the constant upper airway resistance and RA,i is the bronchial inlet resistances of the alveolar compartment .  is the total number of alveolar compartments (for the results in this paper,   = 100).The total flow of air entering the SD at time t k is calculated by During the inhaling phase,   ≥ 0, while in the exhaling phase   < 0. During gas movement in the SD, the fractions of gases in the layer  of the SD, F , ( = 1, … ,   ) is updated based on the composition of the total flow,   , and the current composition of F , .If   ≥ 0, then air starts filling from the top layer ( = 1) to the bottom layer ( =   ); and vice versa for   < 0.
The volume of gas x, in the  ℎ alveolar compartment (v , x ), is given by: In (3), x is any of the five gases (O2, N2, CO2, H2O or α).The total volume of the  ℎ alveolar compartment,   is the sum of the volume of the five gases in the compartment.
For the alveolar compartments, the tension at the centre of the alveolus and at the alveolar capillary border is assumed to be equal.The respiratory system has an intrinsic response to low oxygen levels in blood which is to restrict the blood flow in the pulmonary blood vessels, known as hypoxic pulmonary vasoconstriction (HPV).This is modelled as a mathematical function, resembling the stimulus response curve suggested by Marshall [11], and is incorporated into the simulator to gradually constrict the blood vessels as a response to low alveolar oxygen tension.The atmospheric pressure is fixed at 101.3 kPa and the body temperature is fixed at 37.2°C.
At each t k , equilibration between an alveolar compartment and the corresponding capillary compartment is achieved iteratively by moving small volumes of each gas between the compartments until the partial pressures of these gases differ by <1% across the alveolar-capillary boundary.The process includes the nonlinear movement of O2 and CO2 across the alveolar capillary membrane during equilibration.
In blood, the total O2 content (CO2) is carried in two forms, as a solution and as oxyhaemoglobin (saturated haemoglobin): In this equation, SO2 is the haemoglobin saturation, Huf is the Hufner constant, Hb is the haemoglobin content and O 2 is the O2 solubility constant.The following pressure-saturation relation, as suggested by [12] to describe the O2 dissociation curve, is used in this model: S O2 is the saturation of the haemoglobin in blood and P O2 is the partial pressure of oxygen in the blood.
As suggested by [13], P O2 has been determined with appropriate correction factors in base excess BE, temperature T and pH (7.5005168 = pressure conversion factor from kPa to mmHg): The CO2 content of the blood (CCO2) is deduced from the plasma CO2 content (CCO2plasma) [14] by the following equation: where S O2 is the O2 saturation, Hb is the haemoglobin concentration and pH is the blood pH level.The coefficients were determined as a standardized solution to the McHardy version of Visser's equation [15], by iteratively finding the best fit values to a given set of clinical data.The value of C CO2plasma is deduced using the Henderson-Hasselbach logarithmic equation for plasma CCO2 [16]: where  2 is the plasma CO2 solubility coefficient and pK' is the apparent pK (acid dissociation constant of the CO2 bicarbonate relationship).P CO2 is the partial pressure of CO2 in plasma and '2.226' refers to the conversion factor from millimoles per liter to ml/100ml.[16] gives the equations for  2 and pK' as: For a given pH, base excess (BE), and haemoglobin content (Hb), HCO3 is calculated using the Van-Slyke equation, as given by [17]: The capillary blood is mixed with arterial blood using the equation below which considers the anatomical shunt (ℎ) with the venous blood content of gas x (C v, x ), the non-shunted blood content from the pulmonary capillaries (C cap, x ), arterial blood content (C a, x ), the arterial volume (v  ) and the cardiac output (CO).
The peripheral tissue model consists of a single tissue compartment, acting between the peripheral capillary and the active tissue (undergoing respiration to produce energy).The consumed O2 (VO2) is removed and the produced CO2 (VCO2) is added to this tissue compartment.As per the alveolar equilibration, peripheral capillary gas partial pressures reach equilibrium with the tissue compartment partial pressures, with respect to the nonlinear movement of O2 and CO2.Metabolic production of acids, other than carbonic acid via CO2 production, is not modelled.After peripheral tissue equilibration of gases, the venous calculations of partial pressures, concentrations and pH calculations are done using comparable equations to those above.
A simple equation of renal compensation for acid base disturbance is incorporated.The base excess (BE) of blood under normal conditions is zero.BE increases by 0.
where   =     2 /200000 and   = 0.2  /  Equation ( 16) determines the alveolar pressure   (as the pressure above atmospheric in cmH2O) for the  th compartment of NA number of alveolar compartments for the given volume of alveolar compartment,   () in millilitres.The alveolar compartments are arranged in parallel and interact with the series dead space with respect to the movement of gases.The use of the square of the difference between   and Vc causes alveolar pressure to increase at volumes below Vc, leading to exhalation and a tendency to "snap shut" (mathematical note: the pressure with respect to volume is thus a U-shaped curve) [18].
(per alveolar unit, in cmH2O) represents the effective net pressure generated by the sum of the effects of factors outside each alveolus that act to distend that alveolus; positive components include the outward pull of the chest wall, and negative effects include the compressive effect of interstitial fluid in the alveolar wall.Incorporating   in the model allows us to replicate the situation of alveolar units that have less structural support or that have interstitial oedema, and thus have a greater tendency to collapse.A negative value of   indicates a scenario where there is compression from outside the alveolus causing collapse.The parameter   is a scalar that determines the intra-alveolar pressure for a given volume (with respect to a constant collapsing volume   ) and is dependent on the parameter .
The units of   are cmH2O ml -2 .  , which represents the pressure generated by the respiratory muscles acting on the lung, is described in the next section.Finally,   is defined as a "constant collapsing volume" at which the alveolus tends to empty (through Laplace effects) and represents a fundamental mechanical property of tissue and surfactant [18].  is the end expiratory volume of the lungs.
The effect of the three parameters on the volume-pressure relationship of the alveolar compartments can be observed in Figure S2.The nominal values for ( , ,   ) have been determined such that at the end of expiration, the alveolar pressure within the compartment is also equal to zero, i.e. at an pre-set value of functional residual capacity.
We consider each of the parameters mentioned above ( , ,   ) to be different yet essential components for representing a diseased lung that affect the volume pressure relationship of the alveolar compartments.For example, for a given volume   , increasing   increases the corresponding alveolar pressure of the alveolar compartment.When compared to another compartment with a lower   , a larger pressure gradient would be needed to drive air into the compartment; thus effectively the compartment will be behaving as a stiffer lung unit.
Decreasing  , decreases the alveolar pressure such that the pressure gradient (especially during exhaling) forces the air out of the alveolar compartment until the volume of the compartment collapses (  = 0 ml).Note that, in effect, the parameters are influencing the resting volume of the compartments (when the alveolar pressure,   , is equal to zero).
In the model, the airway resistance   is determined by the following equation for N parallel compartments: In the above equation, the lower airway resistance,    , is determined by the following: where  , is the bronchial inlet resistance of the  ℎ compartment, which is defined by: where  0 corresponds to the default bronchial inlet resistance of an alveolar compartment. 0 is set to 10 −7 ×   for a healthy lung (the inlet resistance is higher for a model with more compartments as the volume of each compartment decreases).  is a coefficient of the airway resistance, representing a dynamic change in airway resistance and is determined by the equation: where, ℎ is the pressure in the trachea and   is a value between 5 and 50 cmH2O for the  ℎ alveolar compartment.Additionally, a threshold opening pressure (TOP) at low lung volumes needs to be attained for a collapsed alveolar unit to open.Recruitment is a time dependent process, with different airways recruiting at different times, once the threshold pressure has been achieved [19], [20].The equations within the model are solved iteratively as a discretised system.Each iteration represents a physiological time slice of t (10 ms).The time-dependent recruitment phenomenon is achieved in the model by the introduction of a parameter t o .For collapsed compartments, t o is set to   which represents the time it could take for collapsed alveoli to open after a threshold pressure is reached.Once  ℎ ≥ TOP i is satisfied, the counter t o decrements during every iteration, and triggers the opening of the airway (  = 1) as t o ≤ 0. Otherwise   is set to a high value (10 10 ) to represent a collapsed airway.
(the number of alveolar compartments) is fixed and set by the user (i.e. they do not change during a simulation).Therefore, during a simulation,   , chiefly represents the relatively small changes in inlet resistance during tidal ventilation.Furthermore,  0 are also preset and fixed, and do not change during the simulation.The only change in airway resistance which is dynamic is   which is dependent on the volume   at time(t k ).
Finally, the pulmonary vascular resistance PVR is determined by where the resistance for each compartment  , is defined as 0 is the default vascular resistance for the compartment with a value of 160 •   dynes s cm -5 min -1 , and   is the vascular resistance coefficient, used to implement the effect of hypoxic pulmonary vasoconstriction.
The net effect of these components of the simulation is that the defining, clinical features of acute respiratory disease may be observed in the model: alveolar gas-trapping (with intrinsic PEEP), collapsereopening of alveoli (with gradual reabsorption of trapped gas if re-opening does not occur), limitation of expiratory flow etc.

Section 2 -Modelling Spontaneous Breathing
The pressure at the airway opening (i.e.mouth) is relatively constant at atmospheric pressure, while the alveolar pressure value varies during inspiration and expiration.The active contraction of the respiratory muscles during inspiration generates a negative alveolar pressure relative to the atmospheric pressure.
The variable   , which represents the pressure generated by the respiratory muscles acting on the lung, is modelled as a piecewise function as described in [21] and adapted from [22].The function consists of a parabolic profile during the inspiration phase of the respiratory cycle, representing the progressive increase in pressure exerted by the respiratory muscles, followed by an exponential profile during the expiration phase of the respiratory cycle, characterizing the passive relaxation of the muscles.
During a single respiratory cycle,   at time t k , is calculated as The resulting alveolar pressure   depends on   , as well as on the pressures of the gases within the compartment, the stiffness of the compartment (  ), and parameter   which represents extrinsic pressures acting on the compartment.
Section 3 -Modelling Non-Invasive Respiratory Support

-Modelling High Flow Nasal Cannula Therapy
To accurately recreate the effects of high flow nasal cannula (HNFC) therapy, we need to consider the differing effects on pressure dynamics during both the inspiratory and expiratory phases.When HFNC is applied, the pressure at the nares increases, primarily attributed to the HFNC-induced augmented resistance during exhalation, resulting in PEEP.
To recreate this phenomenon in our model, we calculate the corresponding increase in pressure by multiplying the HFNC flow by the subject's total airway resistance (Raw) (Eq.22).The modelled Raw comprises two main respiratory airway resistances which are located in series.The first resistance, the upper airways (RU, aw), represents the nasal cavity, oral cavity, and trachea whilst the second resistance, the lower airways (RL, aw) represents the bronchi and the bronchioles, and the inlet resistances for multiple alveoli which make up the lungs, placed in parallel.
To simulate the effect of HFNC during the inspiratory phase, we apply an effective pressure at the end of the upper airways (i.e., at the end of the trachea) by calculating the pressure drop which occurs between the mouth and lungs (Eq.23).During the expiratory phase, this flow is reversed, with lung pressure driving the flow and mouth pressure opposing it (Eq.24).

-Modelling Non-Invasive Ventilation
During Non-Invasive Ventilation (NIV), two pressures are provided to the patient in addition to their spontaneous respiratory efforts: positive end expiratory pressure (PEEP) during exhalation, and an additional inspiratory pressure support (PS) is provided above PEEP during inhalation.To accurately simulate the effect of PEEP and PS, the tracheal pressure applied factors in the pressure drop which occurs due to resistances within the airways.
To replicate the effect of breath triggering in the support being received by the patient, PS is only provided to the patient once their lung pressure falls below the set PEEP due to the negative pressure generated by the respiratory muscles.This pressure then lasts throughout the rest of the active inspiratory phase, with the supporting pressure reducing to PEEP at the beginning of the patient expiratory phase.There is assumed to be no asynchrony between the patient's spontaneous respiratory efforts and the support they are receiving.

Section 4 -Calculating Lung Injury Indices
The following subsections list the methodology for the calculation of each of the key indices of lung injury presented in this study.

-Compliance
The respiratory system compliance (  ) is calculated using: where  is the tidal volume, and Δ  is the pressure swing between peak lung pressure and minimum lung pressure.Note that this equation gives the respiratory system compliance, and not the lung compliance, as we are applying   and   directly to the alveolar compartments.
To calculate the lung compliance (  ) we use the elastance equation: where   is the chest wall compliance, estimated to be 4% of the vital capacity [23].
Vital capacity is estimated using the following equations [24]:

-Transpulmonary Pressure
The transpulmonary pressure (  ) was calculated as: where   is the dynamic lung compliance, and Δ is the difference between the current lung volume and the unstressed lung volume (the residual volume,   ).
The residual volume was estimated for each patient using the following equations [25]:

-Pleural Pressure
The pleural pressure (  ) was derived from the transpulmonary and lung pressures (  ): =   −

-Baseline End Expiratory Lung Volume
To calculate the total strain being applied to the lungs, we first need to determine a baseline endexpiratory lung volume against which we will compare the change in inspired lung volume.To estimate this, each of the digital twins had their respiratory rate set to 14 bpm and their muscle pressure adjusted to achieve a tidal volume set to 7 ml/kg of predicted body (to emulate quiet breathing).The end-expiratory lung volume () was then recorded and used as   .

-Total Lung Strain
The total lung strain was calculated as the sum of the static and dynamic strains across the lung [26].
The dynamic and static strains are calculated as follows: =   +

-Total Lung Stress
The total stress applied to the lungs is, by definition, the stress (pressure) applied to the lungs attributed to the difference between the unstressed volume of the lung,   , and the current volume of the lung.Therefore, the total lung stress is equal to the transpulmonary pressure swing (Δ  ) [27].

-Mechanical Power
The energy applied to the lungs per breath is equal to the area enclosed by the pressure-volume curve.
To determine the energy applied per minute (the power), we multiply this area by the respiratory rate and apply a unit conversion factor, as described in [28]:  = (0.5 ×  ×   ) ×  × 0.0000980665

-Driving Pressure
The concept of driving pressure (Δ), whilst being well defined in invasive mechanical ventilation as the difference between plateau pressure (  ) and , is a more challenging concept in spontaneous ventilation due to the absence of a satisfactory   .It has been proposed that Δ can be approximated in spontaneously breathing patients receiving non-invasive respiratory support by using an inspiratory hold to get an estimate of   [29], but this is still an approximation of the true Δ which is generating the inspiratory flow.
Within the digital twins, we are able to calculate the change in lung pressure which occurs during the inspiratory phase as a result of the support and the patient's inspiratory efforts, allowing an accurate calculation of the pressure which is distending the lung.
In the absence of an elevated inspiratory PS (e.g., during HFNC, CPAP, or unassisted spontaneous ventilation), Δ can be calculated as the difference between the end-expiratory lung pressure () and the minimum lung pressure ( , ) (Figures S3a and S3b).In support modalities which deliver an inspiratory PS above the level of PEEP set, the elevated inspiratory mouth pressure must be factored in.To do so, we measure Δ as the difference between the mouth pressure ( ℎ ) and  , (Figure S3c).This is summarised below: The consequence of the above is that when an inspiratory  is present, the support increases the total distending pressure being applied to the lung.This then facilitates a patient reducing their spontaneous respiratory effort whilst maintaining Δ.
While the above two cases are the most likely to occur and hold true in most cases, there is a third case which can occur.If the  provided to a patient above PEEP fails to exceed the patients  due to the presence of intrinsic PEEP (), then the support will fail to augment the Δ.In these cases, most likely to occur when the  provided is very low, the distending pressure of the lung remains the difference between  and  , , and the patient is unable to reduce their spontaneous respiratory effort without seeing a reduction to Δ (Figure S3d).[30], [31], [32], [33] Section 5 -Patient Matching Protocol

-Matching Overview
Each "digital twin" is a computation model which is designed to reflect the measured characteristics of one of the patients from [34].The twin is created by the careful calibration of the model parameters which govern the system functionality by the minimisation of a user defined cost function (see Section 5.3).

-Cost Function
The cost function aims to minimise the error between the simulator output and the clinical observations by adapting the calibration of the model parameters.The data used in the cost function can be found in Table S3.To separate out the effects of the spontaneous respiratory effort and the additional support being provided, the digital twins had their respiratory support removed while their levels of respiratory effort were maintained constant.Reported below is the contribution of the spontaneous respiratory effort to the lung injury indices.

Section 1 -
Figure S1 -Diagrammatic representation of the pulmonary model.

Figure S2 -
Figure S2 -The effect of varying the parameters of Equation (18) on the pressure volume relationship of the model (  = 0) and under mechanical ventilation.

Figure S3 -
Figure S3 -Four driving pressure instances which can occur.The lung pressure is in blue, with the PEEP in orange (dashed) and the PS + PEEP in yellow (dashed).In each case, the endexpiratory lung pressure (EELP) is signified with a cross (x), and the minimum lung pressure (Plung,min) with a plus (+).Each plot has an arrow showing the driving pressure measured.In cases (a) and (b), no PS+PEEP line is present as PS = 0 cmH2O.

Figure S4 -
Figure S4 -(a) The PEEP provided to the success (blue) and failure (red) groups during NIV.(b-c) The respiratory rate and tidal volume recorded from patients used to titrate the level of inspiratory pressure support (PS) provided.The black dashed line indicates the targets.All three panels present he median(s), interquartile range(s), and full range(s).

Figure S5 -
Figure S5 -The % of tidal recruitment observed for the success (blue) and failure (red) groups during NIV.Please note that he median for the success case is 1%, which coincides with the upper bound of the interquartile range, hence why it is not easily visible.
1 per time slice if pH falls below 7.36 (to compensate for acidosis) and decreases by 0.1 per time slice if pH rises above 7.4 (alkalosis).

Table S4
21 6.2 -HFNC Matching Results to Patient Data 22 6.3 -NIV Matching Results to Patient Data