Exhalatory dynamic interactions between patients connected to a shared ventilation device

In this work a shared pressure-controlled ventilation device for two patients is considered. By the use of different valves incorporated to the circuit, the device enables the restriction of possible cross contamination and the individualization of tidal volumes, driving pressures, and positive end expiratory pressure PEEP. Possible interactions in the expiratory dynamics of different pairs of patients are evaluated in terms of the characteristic exhalatory times. These characteristic times can not be easily established using simple linear lumped element models. For this purpose, a 1D model using the Hydraulic and Mechanical libraries in Matlab Simulink was developed. In this sense, experiments accompany this study to validate the model and characterize the different valves of the circuit. Our results show that connecting two patients in parallel to a ventilator always resulted in delays of time during the exhalation. The size of this effect depends on different parameters associated with the patients, the circuit and the ventilator. The dynamics of the exhalation of both patients is determined by the ratios between patients exhalatory resistances, compliances, driving pressures and PEEPs. Adverse effects on exhalations became less noticeable when respiratory parameters of both patients were similar, flow resistances of valves added to the circuit were negligible, and when the ventilator exhalatory valve resistance was also negligible. The asymmetries of driving pressures, compliances or resistances exacerbated the possibility of auto-PEEP and the increase in relaxation times became greater in one patient than in the other. In contrast, exhalatory dynamics were less sensitive to the ratio of PEEP imposed to the patients.


Introduction
Different critical situations can locally generate a collapse in the health care system. In such circumstances, there may be a time window where hospitals could need to provide ventilatory support to a number of patients greater than the existing number of ventilators. This situation ventilation devices are connected (see, for instance, [13]). Furthermore, it is not clear whether the different devices proposed are universal or ventilator type-specific. At present, dual ventilation strategies lack a substantial amount of mandatory reliable clinical studies. The possibility of performing research on actual patients in a clinical environment is not feasible given current circumstances. In the past, in vivo tests were performed with healthy animals [14] or volunteers [15] that may not strictly reflect the ventilation dynamics on severe acute respiratory distress conditions. Nevertheless, these kinds of studies may provide a rationale for action during critical ventilator shortages. More recently, a group of researchers [16] reported an observational study of four patients with acute respiratory failure due to COVID-19. Tests were performed during the course of an hour with satisfactory results using a device that did not allow to control individual PEEP of patients.
This paper considers a dual ventilation device named 'ACRA' (acronym for Enhanced Capacity of Mechanical Ventilators, in Spanish) that works by setting the ventilator on the pressure controlled mode and provides individualized inspiratory pressures and PEEP control for two patients.
An important feature of the circuit ACRA is that special care has been dedicated to obtain a device compatible with controllers of different ventilators. One main drawback is that because of the different valves incorporated to the system, the flow resistance in the exhalatory limbs has been noted to increase. Additionally, the use of a single in line PEEP valve produces a fluid dynamic asymmetry in the exhalation limbs of the patients.
Alterations produced in the exhalation dynamics and the possibility of auto-PEEP of patients connected to a shared ventilation device require an analysis that is absent at present. Even for perfectly identical patients, any shared ventilation scheme will bring consequences in expiratory dynamics for both of them. This, is associated with the confluence of both expiratory limbs to a single node that forces the pressure to be the same at that point for all branches. Therefore, there is a coupling of circuits and individual exhalation dynamics are no longer independent. Furthermore, the exhalatory valve resistance is a parameter that may enlarge the expiration duration for both patients, as the flow traversing this restriction is largely increased when compared to individual ventilation. This effect significantly depends on, among other parameters, the ventilator considered.
The coupled dynamics of exhalation indicate that the phenomena should be determined by the values of different ratios of parameters associated with the patients and with the ventilation requirements.
Given the complexity involved, research of the fluid dynamics enlightens possibilities and limitations of dual ventilation systems. A better understanding of the underlying fluid dynamics behaviour of the ventilator-circuit-patient system can be achieved through numerical models. Furthermore, with adequate fluid dynamics numerical models, data bases associated with different clinical scenarios can be easily enlarged.
Previous research concerning computational modelling of dual ventilation systems has proposed to analyse the delivered tidal volumes. In [9], a simplified model has been used to compare dual naive ventilation and dual ventilation systems operating in volume control ventilation mode and pressure control ventilation mode. A linear lumped resistance-compliance network model [17] or an electrical circuit as an analogue to the ventilator-circuit-patient [12,18] have also been proposed by other groups. In addition, in [3] a criteria to select patients to be connected to the naive dual ventilation system has been proposed based on numerical models, the required level of support on the lower PEEP/higher FIO 2 and values of lung compliances and resistances.
None of these works have analysed in detail the exhalation dynamics. The estimation of the exhalation characteristic times with linear lumped element models, as considered in these works, takes a simple form when two identical patients are ventilated under the same conditions and with the same exhalatory circuits. However, these ideal cases are the exception. As we will show, when asymmetries of expiratory limbs or patients' tidal volumes occur, no simple expression of the characteristic times can be derived with these traditional analogous electric models. At the same time, there are some limitations of the linear lumped element models that have been discussed previously by different authors (see, for instance, [19]). Particular drawbacks of these approaches are that they consider linear behaviours of the different circuit parts and that distributed elements may behave differently from lumped elements.
Hence, to adequately evaluate the coupled dynamics of exhalations, the development of a refined model seems necessary. In brief, this work aims to analyse the potential adverse interactions during exhalation of mismatched patients connected in parallel to a shared ventilator device with individual regulation of PIP and PEEP. The study is based on results issued from a validated numerical model that we have developed specifically for this purpose.
The rest of this article is organized as follows. In section 2 we describe the model and its validation procedure. In section 3 we show and discuss results considering an ACRA device and a range of compliances and resistances of patients compatible with Acute Respiratory Distress Syndrome. Then, in section 4 we discuss the limitations of our model. Finally, we summarize our main findings in section 5.

Methods
Experimental and numerical tests were carried out to study the possible interaction effects between patients of a shared pressure-controlled ventilation device. While a numerical model is a convenient tool that can accurately characterize the behaviour of many complex system, its performance must be validated against experimental tests. Hence, we first present the laboratory benchmark that allowed the validation and calibration of the model which is analysed afterwards.

Experimental setup
Dual ventilation experimental data was acquired using an 'ACRA' device that was connected to an intensive care unit (ICU) ventilator in pressure-controlled mode. By using valves the device allows accurate control in inspiratory and expiratory pressures on both patients.
As opposed to a system in which a single patient is connected to a ventilator, the inspiratory pressure and PEEP are not directly read from the ventilator monitor. These values are determined by means of setting parameters in the ventilator monitor and manipulating the adjustable valves of ACRA within each patient's circuit. This enables the user to independently modify patient driving pressures.
In this study we used two ACCU LUNG Precision Test Lungs (Fluke Biomedical) and two SmartLung 2000 (IMT analytics) lungs to mimic patients. Also, measurements of pressures and flow rates were acquired through orifice plate differential pressure sensors and a FluxMed monitor (MBMed).
A simplified diagram of the experimental set up is shown in Fig 1. Further constructive details of ACRA can be found in Appendix A.
Typical experimental output signals of pressures and flowrates are presented in Fig 2. Probes placed close to the ventilator and close to each mechanical lung allow to characterize the main fluid dynamics behaviours of the inspiratory and expiratory cycles.
While keeping the output of the ventilator constant and adjusting the constriction of pinch valves or the in-line PEEP valve, different signal waveforms can be observed. When closing the gap of the adjustable pinch valves, pressure values at the end of the inspiration decrease, meanwhile turning the knob of the in-line PEEP modifies the signals of expiratory pressure. Fig 2 shows that waveforms of lung pressure do not achieve the plateau pressure of the ventilator output. The same Figure also illustrates that the flowrate signals exhibit a decelerating flow waveform pattern with a nonzero value at the end of inspiration. The existence of this nonzero flow value at the end of the inspiration and the associated pressure drop produced by its passage through the valve explains why the pressure does not attain the plateau pressure. Fig 3 illustrates the effect of adjusting the constriction gap of one patient pinch valve when keeping the output of the ventilator constant. We observe that by closing this gap waveforms are modified, and the pressure at the end of the inspiration decreases. The Figure shows that closing the gap not only modifies the inspiration curves, as would be expected, but also the exhalation curves as a consequence of alterations of the delivered tidal volume to the patients. In this work, we will focus our study on these dynamics through the models we next describe.

Numerical model
The numerical model we next describe enables us to undertake a comprehensive parametric study of patient conditions and ventilatory requirements to reveal in which situations exhalation dynamics of patients may mutually interfere, producing auto-PEEP in any of the convalescent persons. This model of the circuit was developed with Matlab Simulink using Hydraulic and Mechanical libraries. The code we developed with this platform is open and available on [https://github.com/ACRA2020]. Simulink is a data flow graphical programming language tool for modelling, simulating and analyzing multi-domain dynamic systems, and the block diagram considered is sketched in the comprehensive scheme presented in Fig 4. The solver used in Simulink was ode23t. The trapezoidal rule is implemented in this solver to integrate the ordinary differential equations. The method features a variable step, and a relative tolerance of 1E-3 was stipulated. It is important to mention that this method is more appropriate than Runge-Kutta methods (as implemented in ode45) when the problem is stiff, which happens to be the case due to the high temporal gradients imposed by the ventilator pressure curves used as inputs. Similar strategies have been used by other authors analysing dual ventilation systems [12,17]. Our scheme uses the circuit of an ACRA system but can be easily adapted to other dual ventilation circuits. It comprises: (a) the ventilator, (b) the common circuit parts: inspiratory, expiratory and short circuit tubing, (c) the individualized inspiratory flow circuit, (d) the individualized expiratory flow circuit, and (e) the patient 1 and 2 mechanical lungs models respectively.
The main features of the model were estimated assuming the experimental values taken from the ACRA description presented in Appendix A. The circuit is composed of: (a) Ventilator Module: Fluid properties and a periodic pressure signal define two hydraulic pressure sources, which account for inspiratory and expiratory phases. The input signal  friction factor f and the square of the flow rate _ V L i .

PLOS ONE
Where, A is the pipe cross-sectional area, D the Pipe hydraulic diameter, L the pipe geometrical length and ρ the fluid density. As f also depends on roughness, f = f(Re, ε/D), we considered ε ' 0.2 mm. The compliance of the tubings was considered assuming elastic flexible pipe walls. Considering that in a circular pipe dV dP ¼ LpD (c) ACRA Module: A bifurcation of the inspiratory circuit of the ventilator determines in the ACRA two different branches. Each branch has a one-way valve to force the adequate flow sense and a pinch valve that controls the flow rate. Pressure losses generated by these elements were computed using localized flow resistances through: K i is the localized loss coefficient for each element, which is assumed to be constant and experimentally determined (see Appendix B). The expiratory circuit is similar to the inspiratory one. The difference between them is that an external PEEP or threshold valve (modelled with a check valve) was used to control the PEEP value for patient 2. The pressure loss generated by the PEEP valve was modelled using Eq (2), where K i = K PEEP is the PEEP valve loss coefficient experimentally determined (see Appendix B).  (e) Patients Module: A mechanical system that consists of a spring, a damper and a piston was used to represent the lungs of patients 1 and 2 respectively. As they are often studied under electrical model analogies, the equivalence of both approaches is illustrated in Fig  6. Let us remark that we introduced nonlinear laws for the spring and damper elements, attending the particular mechanical lungs properties we used in the experiments (see Appendix B). Position and velocity of the piston (mechanical parameters) correspond to volume displacement and flow rate on the hydraulic circuit.
In Appendix B we include further details on how the different model parameters are determined.
In Fig 7 we show a validation of the model with typical results for tests performed for two different driving pressures, 20 and 10 cm H 2 O. We confirm a remarkable agreement between the data obtained with the model and measurements.

Results and discussion
The characteristic exhalation times of each patient (τ i ) are used to analyse how the volume of gas in the lungs evolve with time when the fluid is released. When exhalation follows an exponential decay with time, 4τ i corresponds to the amount of time needed to empty 98% of the volume. Considering that the available time for the expiratory cycle is t ex = 1/BR − t ins (with breathing rate BR expressed in min −1 ) and that to empty the lungs of patient i the time required is '4τ i , the result is that to avoid auto-PEEP 4τ i � (1/BR − t insp ). In other words, if the ratio 4τ i /t ex is greater than 1, the lungs will not be fully empty at the end of expiration and auto-PEEP phenomenon may occur. For the sake of clarity, we will use this nondimensional ratio defined as t � i ¼ 4t i =t ex to analyse the characteristic times of the patients during the expiration cycle.
To analyse exhalatory dynamics, it is typical to consider an electric analogy with linear lumped elements in which volume is associated with electric charge, flow rate to electric current, pressure to voltage, friction losses to electrical resistance, and compliance to capacitance. In Appendix C a simplified linear lumped element model is presented. In this case, some basic assumptions such as neglecting pressure losses generated on the one-way valves, on the PEEP valve and on the tubing are made. Also, tube compliances effects are not considered in this simplified model.
A refined model including a linearization of resistances and compliances of the missing components of the circuit, does not necessarily improve the accuracy in predicting time scales. This is due to the high nonlinearities of the components and to the large range of volumes and flows involved in the analysed phenomena.
Taking into account these considerations, the simplified analytical model proposed should be considered as a simple tool that provides reference values and that helps to unveil some relevant aspects of the exhalation dynamics in some specific situations.
For instance, the characteristic exhalation time can be easily estimated when both patients' respiratory parameters, circuit, and driving pressures are identical and the exhalatory valve flow resistance is negligible. In this case, the expiratory dynamics of both patients are coincident and the characteristic time can be determined with the conventional expression τ A = RC where we denoted with R the total resistance of patient and expiratory circuit, and C the total compliance (sub-index 'A' refers to analytical solution).
When incorporating the effect of the exhalatory valve flow resistance in the analysis for these patients, it is easy to demonstrate that τ A = (R + 2R V )C. Note, that this time is larger than the one where the patient is treated individually with the same ventilator (τ A,I = (R + R V )C) (sub-index 'I' refers to individual ventilation). In consequence, connecting identical patients in parallel will always increase the time required for patients to empty the lungs, even in the quite favourable situation of identical patients. In general, when asymmetric situations take place, only a lower limit can be obtained for each patient considering uncoupled dynamics with τ A,i � R i C i for patient i.
A simple nondimensionalization of the simplified model equations, presented on Appendix C, indicates that the nondimensional characteristic times of each patient can be expressed in terms of: where R i and C i represent the airway resistance and lung compliance of patient i (assumed constant and independent of flows and volumes), R V the ventilator resistance, PEEP i the PEEP value of patient i, and DP i the exhalatory driving pressure defined as DP i = P i − PEEP i (where P i is the pressure of patient i at the initiation of exhalation). This pressure difference is intimately related to the traditional driving pressure, more popularly known in the context of mechanical ventilation inspiratory dynamics. The limitations of the analytical solutions do not allow to accurately determine the times required for exhalation. For a given ventilator and circuit, one possibility to determine the functions of Eq (3) could be to rely on an experimental approach. However, the number of tests to be performed would be very large.
Therefore, in this context, the use of numerical modelling appears as an interesting and versatile tool to perform the parametric analysis we next describe.

Numerical estimation of characteristic times
In this section, we propose to numerically evaluate the incidence of the different ratios of Eq (3) on the characteristic exhalation time of each patient. For this purpose, we apply the discussed numerical model, and for each test, the ventilator parameters were set to the typical values presented in Table 1.
In all calculations, the driving pressure of patient 1 was set to a fixed value of DP 1 = 15 cm H 2 O and resistance and compliance for patient 1 were respectively R 1 = 5 cm H 2 O (L/s) −1 and  exponential decay function VT i ðtÞ ¼ a þ be À t=t N;i where sub-index 'i' = 1,2 indicates patient number and a and b are the fitting parameters. The results of the linear lumped element model suggests that in real systems, the exhalation should be expressed with a series of exponential functions (see Eq (8)). However, the results of the numerical model shows that this refinement seems unnecessary to capture the essential physics of the problem. With this model for all the studied cases, the fitting of the results with a single exponential function gave a coefficient of determination R 2 higher than 0.99. Hence, we consider that the use of a single time constant for each patient is representative enough to describe the involved exhalatory dynamics.
In the next sections, a parametric study of the characteristic times for the different ratios of parameters will be discussed. To this end, we will consider the nondimensional ratio calculated using either the numerical model t � N;i ¼ 4t N;i =t ex or the simplified linear lumped element model t � A;i ¼ 4t A;i =t ex . Identical patients with identical ventilation conditions. In this section two identical patients with the same conditions of ventilation will be analysed. ventilator exhalatory resistance results in an increase of the characteristic times for both patients. We can observe that even though patients have the same compliance, resistance, ventilation conditions, and identical circuit check valves, their characteristic times are different. This is because the presence of the PEEP valve on one of the patients' expiratory branch. This valve generates an asymmetry of the circuits and therefore differences on the characteristic time ratios. The associated valve loss coefficient in the system considered is K PEEP = 0.0012 cm H 2 O (L/min) −2 . This in-line valve introduces a pressure loss in the circuit of patient 2 even when adjusted to 0 cm H 2 O and therefore not altering the PEEP value established from the ventilator. It is of interest to compare these results with the characteristic times when a single patient is individually connected to the ventilator (estimated with τ A,I = (R + R V )C) and with the one provided by the linear lumped element model (estimated with τ A = (R + 2R V )C).
As we can see in Fig 9, the characteristic time ratios t � N;1 and t � N;2 are much higher than the characteristic time when only one patient is connected to the ventilator. The reason for this result is that when two patients are connected to the ventilator, the pressure loss on the ventilator exhalatory valve increases drastically due to the rise of the total flow (now . However, this effect may be largely mitigated when exhalatory valves of the ventilator exhibit lower values of flow resistance.  When the PEEP valve resistance is negligible, the symmetry is complete and t � N;1 ¼ t � N;2 . However, as the PEEP valve resistance increases, patient 2 shows an increase in t � N;2 or, in other words, there is more restriction on the patient's ability to expel air from the lungs. In contrast, the effect on patient 1 characteristic time is negligible (its value remains almost constant for the simulated cases).
General case. In Fig 10, the characteristic time ratios t � N;1 and t � N;2 for different PEEP ratios (columns) and driving pressure ratios (rows) are presented using a colour map for a R V /R 1 ratio equal to 0.3. A dotted line is also depicted in each Figure separating the auto-PEEP region (characteristic time ratio greater than 1) and the non auto-PEEP region. Also, the same t �

PLOS ONE
depending on the ratio of driving pressure and the asymmetry of patients, auto-PEEP may take place in one or in both patients. The patient with the in-line PEEP valve is more prone to auto-PEEP than the other. The auto-PEEP of any patient is not quite sensitive to the ratio of PEEP considered. The risk of auto-PEEP for both patients is higher when the compliances of the second patient are higher. The risk of auto-PEEP also increases when the ratio R V /R 1 becomes greater for every combination of patients and ventilation parameters.
To characterize the changes produced on the characteristic times due to asymmetries of the delivered driving pressure, it is convenient to consider the coefficient t � the case corresponding to equal driving pressures). Fig 11 depicts the values attained by this coefficient for different resistance, compliance and driving pressure ratios for R V /R 1 = 0.3. We included in the second row the case of equal driving pressures as a reference (in this case, the values of the coefficient are constant and equal to 0). The first row of the Fig 11 shows the case of a 25% decrease on driving pressure 2 and the third row the influence of a 50% increase on driving pressure 2. In all these cases it can be observed that the asymmetries of driving pressure significantly alter the characteristic time of both patients. The effect is more pronounced for the patient with the lower driving pressure. These results can be explained considering that when driving pressure is increased, tidal volume also increases. Consequently, the combined effect of pressure losses and the coupling associated with flows gathering at the node result in more difficulties to empty the lungs of the patients.

Relevance of this study to clinical guidelines
Regulation Agencies and Medical associations in different countries have stated that alternative ventilation devices proposing pressure control modes of ventilation should be able to control PIP and PEEP [12,21]. Taking into account these requirements, strategies of dual ventilation based on the use of simple splitters (naive sharing circuits devices) have been legitimately questioned because mismatch of resistances and compliances may lead to inadequate ventilation of patients. However, previous works have shown that incorporating to the naive sharing circuit devices variable resistances, it is possible to individualize the tidal volume of each patient through controlling the specific PIP delivered to each patient. Furthermore, by incorporating in-line PEEP valves, it could also be possible to individualize the control of PEEP delivered to each patient. However, this last control may also depend on an adequate matching of patients because of possible auto-PEEP. Our study shows that the risk of auto-PEEP is quite limited for identical patients. For similar patients, to mitigate the possibility of auto-PEEP, it is beneficial to connect the patient with lower compliance to the branch that incorporates the in-line PEEP valve. The larger the asymmetries in patient compliances and driving pressures, the larger the probability of auto-PEEP. Moreover, our study shows that the exhalation characteristic times of the patients is not quite sensitive to differences on the ratios of PEEP values.

Limitations of this study
The numerical model we developed here allows the user to consider nonlinearities in compliances and resistances and shows great agreement with in vitro tests. However, in our study, some of these parameters like the values of exhalation resistances and compliances of each patient have been considered as constant with the flow. This may be quite simplistic. Hydraulic resistance for instance in the endotracheal tube will be constant only for laminar flow, a regime that usually takes place only at the end of the expiratory cycle. Also, the physiological relationship between pressure and flow in patients is often nonlinear (see [22,23]). In addition, lung compliance is nonlinearly dependent on ventilation pressure due to volume-dependent tissue stiffness [24].
Taking into account previous research [25] we also considered that the resistance of the exhalatory valve was linear with the flow. Most currently available ventilators have a servocontrolled exhalation valve with a variable opening controlled by the ventilator. Usually, for a fixed position of the valve opening, the pressure drop through the valve depends quadratically with the flow DP ¼ K v _ V 2 . When the valve opening changes with time, the law that links the instantaneous pressure drop with the flow can become complex over time. Particularly, the closure dynamics of the exhalatory valve may differ with the type of ventilator considered. Therefore, assuming a universal linear link between instantaneous flow and instantaneous pressure drop may also be simplistic.
The analysis we have performed considers the characteristics of the pinch valves used in the ACRA device and the characteristics of the exhalatory valve of a Nellcor Puritan Bennet 760 Ventilator. Even if the results here obtained are modified for use with other ventilators, the observed trends should be the same. However, results of this parametric study may not reflect all situations found in clinical care.

Conclusions
In this work, we analyse for the first time the consequences in exhalatory dynamics when two patients are connected to a dual ventilation system. A simplified linear lumped element model was used to obtain reference values and to identify the set of relevant nondimensional parameters related to patients and ventilation requirements involved in the study.
We performed a parametric analysis with a computational model, distinct in its ability to include nonlinear behaviours and distributed parameters. The remarkable agreement between results issued from the model and our experimental results with test lungs enabled us to validate the model. However, the 1D model developed required an adequate evaluation of the valves' flow resistances (circuit and ventilator) to produce reliable results. Our study shows that dual ventilation devices produce larger patient exhalation settling times than they would have with individualized ventilation. As a consequence, patients connected to shared ventilation may experience auto-PEEP. However, this effect can be mitigated with the adequate matching of patients. The presented results correspond to a specific shared ventilation device (ACRA) used with an specific ventilator (Nellcor Puritan Bennet 760). However, the developed open source code of this article, provides the opportunity to extend the analysis of auto-PEEP risk when connecting patients to other ventilators and shared ventilation devices with simple parameter characterizations. At present, our results are based on a computational model that has not been validated against any preclinical or clinical data. Hence, clinicians should not look to our study as an exact estimate of patients exhalatory dynamics.
Appendix A: ACRA dual ventilation system description Fig 12 shows 3D renders of the ACRA device. This apparatus has been tested in two paired breathing simulators (ASL 5000, InGMAR Medical, Pittsburgh, PA) and also with two paired pig models whose lungs were manipulated to mimic severe acute respiratory distress conditions. We performed a 2-hit experimental ARDS model. Repeated lung lavages with normal saline (30 mL kg −1 at 37 degree Celsius) were performed until PaO 2 /FIO 2 = 200 mm Hg at PEEP 10 cm H 2 O. The animals were then submitted to 120 min of injurious mechanical ventilation at PEEP = 0 cm H 2 O, VT = 15 mL kg −1 , breathing rate BR = 12 min −1 , an inspiratoryexpiratory ratio of I:E = 1/2, and a FIO 2 = 1.0. Thereafter, baseline ventilatory settings were restored; the success of the model was confirmed with a PaO 2 /FIO 2 = 200 mm Hg, and lung ultrasound imaging showing bilateral atelectasis in the dependent lungs.
The ACRA system incorporates 700 cm 3 including the volume of the tubing from ventilator to the interface (300 cm 3 = 2 × 0.4 m long 22 mm flextube). The two breathing circuits (one for each patient) incorporate 2400 cm 3 (2 × 1.6 m long 22 mm flextube breathing systems). The total volume of tubings and the ACRA system is 3100 cm 3 .
Each patient's inspiratory pressure is established by introducing a controlled pressure drop from the ventilator output P vent . This pressure reduction occurs mainly in the adjustable flow restrictor (pinch valve) located on each patient's inspiratory circuit.
Changing the value of the gap constriction modifies the flow rate/pressure drop ratio. The gap value in this constriction is adjusted manually by turning a valve knob.
The peak inspiratory pressure (PIP) in the patient is determined by the value of the flowrate at the end of the inspiratory cycle, and it will be lower than the pressure imposed by the ventilator. Therefore, the PIP is not only determined by P vent and the flow restriction but also by the duration of the inspiratory cycle.
Two analog manometers are incorporated in the device to determine the pressure values downstream of each pinch valve. These devices allow to measure each patient's driving pressures. When air is still flowing at the end of inspiration, a ventilator inspiratory pause will be required to read plateau pressure.
PEEP value is established by the ventilator for the first patient. Conversely, the second patient circuit has an adjustable in-line threshold valve (PEEP valve). The threshold value established with this valve in its expiratory circuit is the sum of the PEEP value shown in the monitor and the one set in the adjustable valve. The reading of the threshold value established by the position of the PEEP valve knob is not required. The ACRA manometers directly allow the user to read the values of PEEP, and the position of the knob can be adjusted from this reading.
The set of one-way valves prevent the flow from moving in unwanted directions. The oneway and threshold valves of the expiratory circuit should be selected to produce low values of pressure drop and seal adequately at low and high pressure differences.
A short circuit tubing between inspiratory and expiratory limbs, also proposed by [13], ensures compatibility of the system with different ventilator controllers. As the ACRA system introduces pressure drops in the respiratory circuit, unknown to the ventilator controller, the feedback loop that enforces the different goal pressures throughout a respiration cycle might spiral out of control if the system never reaches its target. This is the case while measuring the target PIP from the expiration circuit, as the pressure drop introduced by the ACRA imposes an effective PIP which might be below target by a noticeable margin. Depending on the monitoring enforced, this could either raise an alarm, or simply lead to a feedback loop in the inlet regulation providing a higher flow rate and, eventually, delivering excessively high volumes to the patients. Also, a high discrepancy between inspiratory and expiratory pressures could lead to monitoring errors or false alarms [13].
Therefore, to establish an adequate protocol independent of both the ventilator considered and the feedback pressure probe location, it is beneficial to use short-circuit inspiratory and expiratory circuits as indicated in Fig 1(green line in single-line diagram) with a small size tubing (i.e. 6 mm in diameter). With this scheme, the ventilator feedback control system senses a single 'virtual' patient with a small resistance (mainly that of the short-circuit tubing) and a large compliance (mainly those of both patient lungs). Hence, the controller will target the parameters fixed on the monitor and adjust internal variables to obtain the predetermined goal of inspiration as a consequence of any disturbance with this patient configuration.
In [18] it has been emphasized the need to monitor the delivered tidal volume to each patient. The reading on the ventilator screen of the tidal volume with the ACRA corresponds to the value issued from the addition of volumes of both patients. To determine individual values, it is always possible to interpose in the inspiratory branch of a patient a flow sensor. Also, individual tidal volumes can be directly read on the ventilator screen by clamping the endotracheal tubing of one patient during a few cycles.

Appendix B: Model parameters estimation
In our experimental set up, the values of lung resistance and compliance depend nonlinearly on the flow and volume respectively. This means that the compliance of both patients depends on the volume of each lung (C i = f(V i )) and the resistance of both patients depends on the flow . These functions were incorporated into the model extracted from the SmartLung 2000 characteristic curves presented by the manufacturer in the datasheet. A set of system parameters were obtained by fitting the results of laboratory experiments carried out for different flow conditions. The pressure loss curves of filters, one-way valves, and PEEP valve used in the model are shown together with the experimental measures in Fig  13, and the pressure loss coefficients are presented in Table 2.
In all these fittings a quadratic dependence of pressure loss with flowrate was considered. It is not easy to directly determine the expiratory valve characteristics of the ventilator. The governing law of this valve can take different complex forms and in many cases the geometry must be taken into account to correctly characterize it [26][27][28]. In our model, a linear pressure flow relation was considered (DP ¼ R V _ V ). Constant R V was determined through a minimization between pressure experimental measurements P exp aw i and the corresponding numerical output of the signal P num aw i ðR V Þ. Experimental data was obtained for each patient circuit and for the expiratory circuit to ventilator. Flow rates and pressure signals were measured through orifice  plates and differential pressure sensors connected to a FluxMed monitor. Numerical outputs were produced considering different values of (R V ). To determine the more appropriate value of this parameter, we defined a functional J, where P exp aw 1 , P exp aw 2 and P exp esp are patient 1 experimental airway pressure, patient 2 experimental airway pressure and experimental expiratory pressure respectively. Also, P num aw 1 , P num aw 2 and P num esp are patient 1 numerical airway pressure, patient 2 numerical airway pressure, and numerical expiratory pressure respectively all as a function of R V . Functional J is plotted against R V in Fig 14 for  approach. In this simplified linear model, resistances and compliances are considered constant in time and independent of the involved flows or air volumes.
Kirchhoff circuit laws lead us to a system of differential equations presented in Eqs (5) and (6), where we identify flow rate of the patients with V i (with i = 1, 2) and dots indicate time derivative. Denoting PEEP 2 = ΔPEEP 2 + PEEP 1 and expressing the system of equation in a matrix form we get, This equation system � � C: � V þ � � R: _ � V ¼ � PEEP is a system of first order linear nonhomogeneous differential equations. Therefore, the general solution of the problem will be the sum of the homogeneous and a particular solution.
This solution takes the form, where the scalar λ i and vector � A i are the eigenvalues and eigenvectors of matrix ( � � C À 1 � � � R) respectively. The eigenvalues λ i can be computed as l i ¼ 1 2 ðb � ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi b 2 À 4c p ), where b = C 1 (R 1 + R V ) + C 2 (R 2 + R V ) and c ¼ C 1 C 2 ðR 1 þ R V ÞðR 2 þ R V Þ À C 1 C 2 R 2 V . The general expression for the eigenvectors can be seen on https://github.com/ACRA2020/Analytical-Solution. Finally, scalars α and β depend on the values of tidal volumes at the end of inspiration.
In order to obtain a tractable analytical solution, the assumption of null flow at the end of inspiration is made. Then, considering pressures at the end of inspiration P i as the initial conditions of expiration, the constants α and β can be determined. It is worth mentioning, that in this case, the exhalatory driving pressure is equal to the driving pressure. For a more general case, the exhalatory driving pressure DP i for each patient is the difference DP i ¼ P i À PEEP i (with P i the pressure of patient i at the initiation of exhalation) and, Note that this maximum pressure difference is, in this case, intimately related with the driving pressure defined as the plateau airway pressure minus PEEP.
These last equations indicate that constants α and β depend on the driving pressures and compliances of each patient. Therefore, this analysis shows that the speed at which lungs go from an initial state (that is, function of the pressure at the end of inspiration) to the final state depends on R 1 , R 2 , R V , C 1 , C 2 , DP 1 and, DP 2 .
With the aim of identifying the relevant physical quantities that dictate the exhalatory fluid dynamics, we nondimensionalize the simplified model equations with patient 1 parameters.
We define the nondimensional volume as, v 1 = V 1 /(C 1 DP 1 ) and v 2 = V 2 /(C 1 DP 1 ) and the nondimensional flow as, Taking into account these definitions, the homogeneous system of Eqs (5) and (6) can be written in a new and equivalent system of first order nondimensional linear homogeneous differential equations as: By simple inspection of Eqs (10) and (11), the homogeneous system of equations presented is governed by three groups of parametric relations, R 2 /R 1 , C 2 /C 1 , and R V /R 1 The same reasoning can be applied for Eq (9). Then, it can be concluded that the general solution of Eqs (5) and (6) depend on the following ratios, R 2 /R 1 , C 2 /C 1 , R V /R 1 , DP 2 /DP 1 , PEEP 2 /PEEP 1 .
The use of patient 1 parameters to undertake the nondimensionalization of the equations establishes that the expiratory characteristic times obtained are referred to the time constant R 1 C 1 . However, to facilitate the analysis without loosing generality, we preferred in the manuscript to present the expiratory characteristic times referred to the available exhalatory time t ex .