Simulink Modelling For Simulating Intensive Care Mechanical Ventilators

This paper proposes a modelling approach for simulating mechanical ventilators for intensive care units (ICUs). The shortage of ventilators during the coronavirus disease 2019 (COVID-19) pandemic has focused attention on their design and performance. The proposed modelling approach consists in using the Mathworks®Simulink software tool and the SimScape Fluids (gas) library, so as to use well-established subroutines to simulate all the pneumatic components of typical ventilators for ICUs, such as the pressure reducing valves, pressure relief valves, check valves, tanks, ON\OFF and proportional directional valves, etc. The patient is simulated by setting the values of lung compliance and pressure losses occurring in the trachea. The proposed modelling approach is used in this paper to simulate a pneumatic scheme employed in some commercial ventilators. The model allows a very accurate prediction of fundamental parameters, such as the inspiratory flow rate, the inspiratory pressure, the end-expiratory pressure. Since the software interface is user-friendly, it can easily be used by manufacturers to correctly choose the geometrical and operating parameters of the components during the design stage or to assess different scenarios.


Introduction
A mechanical ventilator is a machine that provides assisted ventilation to patients that are totally or partially unable to breathe on their own because of critical conditions due to illness or surgical procedures [1]. Mechanical ventilators are fundamental elements for intensive care units (ICUs); the importance of ICUs has recently been elevated due to the worldwide pandemic of coronavirus disease 2019 (COVID- 19), since they are necessary to sustain critically ill patients who are not able to breathe naturally after becoming infected by SARS-CoV-2 [2].
An ICU mechanical ventilator is a complex machine, which can work with different modes of operation and needs several parameters to be set, such as the respiratory rate (RR), pneumatic scheme using the International Organization for Standardization (ISO) symbols. Mechanical ventilators for intense care units are usually connected to wall outlet supplies of compressed air and oxygen, which are prone to significant pressure variations ranging from 2 to 7 bar gauge [1]. The air and oxygen enter the device via the gas inlet connections (0). In the first part of the ventilator, the air and the oxygen follow parallel paths. After flowing through a filter (1), used to remove impurities, and a check valve (2), needed to avoid reverse flow, both the air and the oxygen enter a pressure reducing valve (3). The role of the two pressure reducing valves is to decrease the pressure of the two gases coming from the wall outlet supplies; the desired value of the pressure downstream of the two pressure reducing valves, measured by the pressure sensors (4) and being usually of the order of 500 mbar gauge, can be adjusted manually or using a control signal. The blend of the air and oxygen is obtained through the air metering valve (5) and the oxygen metering valve (6), which are usually directional proportional valves [27]. The opening degree of these valves can precisely be adjusted by acting on the control signal sent to the solenoids, thus obtaining the wanted mixture of gas with the selected percentage of oxygen. Depending on the application, using high-pressure proportional valves can eliminate the need for the pressure-reducing valves [27].
The mixed gas enters a tank (7), which is needed both for making the mixture of air and oxygen uniform and for smoothing the pressure variation at this point of the circuit. The model Galileo [9] is equipped with a tank having a volume of 5.8 L in which the pressure of the mixture is maintained in the range from 200 to 350 mbar gauge. According to the manufacturer, a further advantage of this large volume is that it is able to supply high peak flows when required [9]. The tank pressure-relief valve (8) has the task of avoiding uncontrolled pressure increase in the tank due to malfunctioning. Therefore, when the pressure in the tank reaches a pre-set maximum value, the valve automatically opens and discharges the gas into the room in which the ventilator is located.
The oxygen concentration in the tank is monitored by an oxygen sensor (9), which generates an output voltage proportional to the oxygen concentration, used to control the opening degree of the inlet metering valves (5) and (6). An alarm can be present to report an unusual concentration of oxygen.
An additional filter (10) can be used to protect the downstream inspiratory valve (11) and the patient against possible contamination particles carried by the gas. In modern ventilators, the inspiratory valve (11), which allows precise control of the gas inspired by the patient, is a proportional valve whose opening degree can be controlled accurately. In some models, such as Galileo [9], this valve is also equipped with a position sensor coupled with a differential pressure sensor measuring the pressure difference across the orifice of the inspiratory valve. In this way, it is possible to calculate the flow of gas through this valve. In other few models (e.g., small transport ventilators and automatic resuscitators), this valve is just a fixed orifice flow resistor that delivers a constant flow. In simple infant ventilators (e.g., Bourns BP-200 and Infant Flow SiPAP device by CareFusion), the inspiratory valve is manually adjusted [1,28].
A pressure relief valve (12) is placed in the circuit in order to avoid uncontrolled pressure increase downstream of the patient, which could harm the patients' lungs. Therefore, when a pre-set maximum value of the pressure is reached (namely, the cracking pressure), this valve automatically opens and discharges the gas into the external environment. A further safety valve (13) is implemented in the circuit: this valve, which is normally open, allows the patient to breath in case of standby or malfunctioning of the system. During normal operation of the ventilator, these valves (12 and 13) remain closed. The pressure at this point of the circuit is measured by the pressure sensor (14).
A humidification device and a heat exchanger (15) can optionally be used to control the humidity and temperature of the gas inspired by the patient.
A flow sensor (16) is positioned close to the patient's mouth in order to measure and monitor the flow of gas inspired and expired by the patient.
The final part of the circuit consists of a pressure sensor (17), an optional additional filter (18), an expiratory valve (19) and a check valve (20). In the Galileo [9] and Dräger Infinity V500 [10] models, the expiratory valve is a proportional directional valve which is used both to adjust the minimum pressure in the patient circuit and to regulate the flow of expiratory gases from the patient to the external environment. The check valve (20) is used to avoid reverse flow in the expiratory circuit.
Some medical treatments need the patient to be provided with some medicine in aerosol form. To accomplish this task, different strategies can be used [1]. In Galileo [9], when this additional function is activated, the directional on/off valve (21) is opened. Afterwards, the compressor (22) increases the pressure of the gas taken from the tank (up to 900 mbar gauge), which can enter a nebulizer jar (23) to blend with the medicine to be supplied to the patient. This method has the advantage of not disrupting the air/oxygen ratio of the gas that Galileo delivers to the patient [9]. Furthermore, during nebulization Galileo automatically reduces the flow of gas through the inspiratory valve, thus compensating for the additional flow of gas delivered to the patient by the nebulizer [9]. Concerning the control system, there can be different control strategies. Among these, one of the most used for critical patients is called "Volume Controlled Ventilation (VCV)", which is a mode of ventilation where a set gas volume (tidal volume) is delivered to the patient [14]. To achieve this target, the inspiratory valve (11) is opened and the flow rate inspired by the patient is measured by the flow sensor (16). When a correct volume of gas (tidal volume) has entered the patient's lungs, the inspiratory valve (11) is closed and the expiratory valve (19) is opened to allow the patient to expire. In this process, the positive end-expiratory pressure (PEEP) is kept higher than a pre-set value to prevent lungs from collapsing.
Another control strategy, which is called "Pressure Controlled Ventilation (PCV)", is a pressure-targeted mode of ventilation, in which the ventilator adjusts the inspiratory flow to keep the maximum airway pressure at the set level during inspiration [14].
In addition, a ventilator should have a spontaneous breathing pressure support mode for those patients breathing to some extent themselves, namely the "Bilevel Positive Airway Pressure (BIPAP)" mode (a non-invasive ventilation mode that provides different levels of pressure when the patient inhales and exhales) or the "Synchronized Intermittent Mandatory Ventilation -Pressure Controlled (SIMV-PC)" mode (a ventilation mode where the patient is allowed to take spontaneous breaths, the machine will assist the patients breathing when a spontaneous breath is taken) [14].
Concerning the operating parameters that a ventilator needs to meet, the Medicines and Healthcare products Regulatory Agency (MHRA) has recently divulged some guidelines and operating parameters to be met [14]. Among these, there are:  The plateau pressure, which is the pressure in the lungs measured at the end of the inspiratory pause, must be less than 34.32 mbar (35 cmH2O) gauge [14].  The mechanical failsafe valve (see pressure relief valve (12)) must open at 78.45 mbar (80 cmH2O) gauge [14].  The Positive end expiratory pressure (PEEP) must be greater than 4.903 mbar (5 cmH2O) gauge, adjustable in at least 4.903 mbar (5 cmH2O) increments.

Simulation code
A simulation code, reproducing the pneumatic scheme of Fig. 2, has been realised in Mathworks ® Simulink using the SimScape Fluids (gas) library [25,26]. The code, whose graphical scheme is shown in Fig. 3, can easily be reproduced and used by manufacturers and is available at the link provided in ref.
[29]. The gas properties are defined in the "Gas Properties" block, allowing the user to select among three gas property models: perfect gas, semiperfect gas, and real gas. The specific gas constant is also set in this block. Only the properties of one gas can be defined; therefore, the gas constant must be chosen to fix the oxygen/air mixture composition (i.e., the percentage of O2 in the gas), which is equivalent to fixing the molar mass of the mixture. The content of humidity in the gas is not simulated.
The air and oxygen sources (0) are simulated using the "Pressure Source" block. This block maintains a constant pressure at its outlet regardless of the mass flow rate through the source.
The "Check Valve" blocks, placed downstream of the gas sources, simulate the check valves (2). With this block, the flow is allowed in one direction only; the cracking pressure, which is the minimum upstream pressure required to open the valve, is set in this block.
Two "Pressure Reducing Valve" blocks are used to simulate the pressure reducing valves (3). The value of downstream pressure that is maintained by this block is defined in the block properties.
Two "Pressure & Temperature Sensor" blocks reproduce the pressure sensors (4) positioned downstream of the pressure reducing valves to monitor the pressure at these points of the circuit. The air and oxygen metering valves (5 and 6) are simulated using the "Variable Orifice ISO 6358" block. This block models a controlled pressure loss from port A to port B based on the ISO 6358 standard, which can represent a valve, an orifice, or a restriction. The orifice opening fraction between 0 (valve closed) and 1 (maximum opening) is set by the control signal at port L. This block allows the valve properties to be set by using the sonic conductance, the flow coefficient or the restriction area of the corresponding valve. The critical pressure ratio of the gas must be defined as well.
The tank [7] is reproduced by the "Constant Volume Chamber" block, by which the tank volume is set. The pressure in the tank is calculated according to the mass, the volume and the temperature in the tank. The sole inlet port is port A, whereas there are two outlet ports: port B, connected to the downstream inspiratory valve (11), and port C, connected to the "Nebuliser".
The inspiratory valve (11) is simulated by the "2-Way Directional Valve" block. The flow rate is calculated according to the ISO 6358 standard. Ports A and B are the gas ports. The valve opening fraction between 0 and 1 is set by the control signal at port S. A positive signal opens the connection between ports A and B proportionally. The valve properties can be set in the block by using the sonic conductance, the flow coefficient or the restriction area of the corresponding valve.
The pressure relief valve (12) is simulated by the "Pressure Relief Valve" block. The valve remains closed when the upstream pressure is less than the maximum pressure set in the block. When the maximum pressure is reached, the valve opens to permit gas flow. The safety valve (13) is reproduced by the "2-Way Directional Valve Block". The valve is closed when the control signal is 0 and opens when the control signal is 1. These two valves remain closed during normal operation. A "Pressure & Temperature Sensor" is positioned at this point (14) to measure the pressure downstream of the inspiratory valve.
A "Pipe" block (15) is implemented to take into account the pressure drops in the circuit due to the presence of the heat exchanger and humidifier, by using an equivalent length and an equivalent diameter. By this block, it is also possible to set the value of the temperature of the gas exiting the heat exchanger.
A "Volumetric Flow Rate Sensor" block (16) is positioned upstream of the patient to measure the flow rate inspired (>0) and expired (<0) by the patient.
The expiratory valve (19) and the check valve (20) are simulated using the already introduced "2-Way Directional Valve" block (the valve opening fraction between 0 and 1 is set by the physical signal at port S) and "Check Valve" block, respectively.
The patient is simulated using the subroutine "Patient" shown in Fig. 4. In this subroutine, the pressure losses occurring in the trachea are reproduced by the "Orifice ISO 6358" block, in which an equivalent sonic conductance, an equivalent flow coefficient or an equivalent restriction area can be defined. The lungs are simulated using the "Translational Mechanical Converter" block (namely, a piston moving inside a cylinder), which is connected in series with the "Spring" block simulating the lung compliance. The piston diameter and the spring stiffness must be chosen in accordance with the lung compliance of the simulated patient. The lung compliance C is defined as follows: where ∆ is the volume variation of the lungs caused by the pressure variation inside the lungs ∆ . Having reproduced the lungs as a cylinder whose piston is counteracted by a spring, the relation between ∆ and ∆ is: where A is the piston area, ∆ is the displacement of the piston and K is the stiffness of the spring. Combining equation (1) and equation (2), one obtains: Equation (3) allows one to calculate K once A has been fixed (or vice versa), in order to obtain the wanted value of compliance. The initial compression of the spring x0 can be set to the match the minimum pressure in the circuit occurring at the beginning of the inspiratory phase, namely, the Positive end expiratory pressure (PEEP): Finally, the dead volume of the cylinder can be calculated as follows: The subroutine reproducing the nebuliser is shown in Fig. 5. The ON/OFF valve (21), the compressor (22) and the nebuliser jar (23) are reproduced, respectively, by the already introduced "2-Way Directional Valve" block (the control signal being equal to either 0 or 1), "Pressure Source" block (which maintains a constant pressure at its outlet) and "Orifice ISO 6358" block (in which an equivalent sonic conductance, an equivalent flow coefficient or an equivalent restriction area can be defined).
Two proportional-integral-derivative (PID) controllers have been implemented in the program: "PID flow" and "PID pressure" (see Fig. 3), respectively allowing the automatic control of the inspiratory flow rate and the automatic control of the positive end-expiratory pressure (VCV mode). These subroutines are not described here for brevity; the reader can find these subroutines at the link provided in ref. [29].
The resulting set of ordinary differential equations, describing the dynamics of the system with respect to time, is solved by the solver Ode 23t with variable time steps (from 0.001 s to 0.1 s) [25,26].

Results
Using the numerical code described in Section 3, the respiratory cycle can be simulated for different control strategies, for different patients' characteristics (i.e., lung compliance and pressure losses in the trachea), and for different geometric characteristics (i.e, size of the valves and of the tank) and operating parameters (i.e., pressures and flow rates) of the ventilator. In this way, several scenarios can be predicted.
Some examples of simulation are now described and discussed. In the simulated control strategy, the VCV mode is used, and the inspiration is triggered by the ventilator itself, with a null patient's effort.
The operating and geometrical parameters chosen for the simulated ventilator are shown in Table 1. The desired percentage of O2 in the gas is taken as equal to 20.9 %, which means that the oxygen metering valve is maintained closed. The nebulizer function is not activated by the operator in this example.
Concerning the geometrical parameters, the sonic conductance at maximum flow rate for the air metering valve, oxygen metering valve, inspiratory valve and expiratory valve is taken as equal to 10 L/s/bar. The opening offset of the inspiratory and expiratory valves is assumed as equal to -10% of the full opening. The tank volume is considered as equal to 6 L, with a temperature inside the tank of 20 °C. The pressure drop caused by the heat exchanger and humidifier is taken into account by an additional equivalent length of 0.5 m and an additional equivalent diameter of 0.01 m. The heat exchanger is assumed to heat the gas up to 27 °C. Table 2, three different patients, having different values of lung compliance, are simulated. The lung compliance is 0.04079 L/mbar (0.04 L/cmH20) for patient 1, 0.06118 L/mbar (0.06 L/cmH20) for patient 2, and 0.05099 L/mbar (0.05 L/cmH20) for patient 3. The pressure losses in the trachea are 5 mbar (5.099 cmH20) at 25 L/min for the three patients. The temperature of the gas inside the lungs is 37 °C for the three patients.

As shown in
The inspiratory time, the inspiratory pause, and the expiratory time are defined by the operator, which means that the opening and closing times of both the inspiratory valve and the expiratory valve are set accordingly by the control system. Patient 1 and patient 2 have the same inspiratory and expiratory time intervals, which differ from those of patient 3.
The target tidal volume (TV) is set by the operator: having set the inspiratory time, the inspiratory flow rate (IFR) becomes the directly controlled variable. The control system (subroutine "PID flow") acts on the opening degree of the inspiratory valve to reach this set point. The operator also sets the wanted positive end expiratory pressure (PEEP) value. The control system (subroutine "PID pressure") acts on the opening degree of the expiratory valve to reach this set point.
The combinations of the target IFR and PEEP values set in the simulations are as follows: IFR =25 L/min (TV=0.5 L) and PEEP = 5 mbar (5.099 cmH20) gauge for patient 1 and for patient 2; IFR =30 L/min (TV=0.5 L) and PEEP = 10 mbar (10.20 cmH20) gauge for patient 3.  Table 3 reveals the parameters entered in the subroutine "patient" to obtain the values of lung compliance and pressure losses in the trachea as well as the initial conditions for the three simulated patients (see Equations 1 to 5).
The values considered in the simulations are consistent with the values present in the literature [15][16][17][18].
The predicted time history of the pressure upstream of the patient (pressure sensor 14) is reported in Fig.6 (a) for patient 1 and in Fig. 7 (a) for patient 2; likewise, the predicted flow rate inspired (>0) and expired (<0) by the patient is shown in Fig. 6 (b) for patient 1 and in Fig. 7 (b) for patient 2.
The comparison between the predicted trends for these two patients, which have equal target values of IFR and PEEP, show that the results are consistent with the expected trends. Indeed, higher peak values are obtained for the pressure upstream of patient 1 and for the expired flow rate of patient 1, due to the lower lung compliance of patient 1. For both patients, the control system has been able to reach, accurately and without overshoots, the set points defined for the inspiratory flow (25 L/min) and for the positive end expiratory pressure (5 mbar gauge).  Table 3. Parameters set in the subroutine "patient" Notably, for all the three patients considered, the initial breathing cycle shows some differences compared to the following breathing cycles. This is due to the initial conditions assumed in the simulations. A better tuning of the initial conditions may lead to a more effective simulation of the initial cycle.  Fig. 9, the predicted trends of the pressure upstream of patient 1 and of the flow rate inspired and expired by patient 1 are enlarged and compared with the corresponding expected trends provided in the literature [30]. The comparison shows that the numerical predictions are in very good agreement with the expected trends.
In addition, the numerical values obtained for the fundamental parameters of the respiratory cycle, such as the peak inspiratory pressure (PIP), the plateau pressure (Pplat), and the positive end expiratory pressure (PEEP), are consistent with the values specified in the guidelines provided by the MHRA [14].  The proposed code is also useful to predict the behaviour of the valves employed in the circuit. In this regard, Fig. 10 shows the behaviour of the inspiratory and expiratory valves during a breathing cycle, with reference to patient 1. It is shown that, as expected, the inspiratory valve is opened during the inspiratory phase and then it is closed during the expiratory phase, ensuring that the inlet gas is not lost during the patient's expiration. However, this valve is opened up to only 14 % of the maximum opening degree, which means that, for the chosen operating conditions of the ventilator, a smaller valve can be employed in order to exploit the entire opening range of the valve. Concerning the expiratory valve, it has been selected correctly in terms of size, since the entire opening range is exploited during the expiratory phase. The maximum opening degree is obtained at the beginning of the expiratory phase, when the pressure upstream of the patient is rapidly reduced from the plateau value down to the PEEP value defined as the set point.

Conclusions
This paper provided a Simulink modelling approach for simulating typical ventilators for ICUs, which is intended to be a tool to be used by manufacturers and start-ups to speed up the design and the production of new mechanical ventilators for ICUs. The proposed modelling approach was used in this paper to simulate a pneumatic scheme employed in some commercial ventilators. The obtained model was thoroughly described in the paper in order to allow one to easily reproduce and use it. The simulations presented in the paper show that the modelling approach is reasonable, providing an accurate prediction of the expected respiratory cycle. The code is also made available at the link provided in ref. [29].
The user-friendly interface, where each component is represented by a block in which the geometrical and operating parameters can be entered, can be used by manufacturers for different purposes, such as:  to assess the effects of different geometrical parameters (such as the sonic conductance, flow coefficient and restriction area of the directional valves, the tank volume, the cracking pressure of the check valves, etc.);  to assess the effects of different operating parameters (such as the supply pressure, the pressure downstream of the pressure reducing valves, etc.);  to assess how the ventilator responds to different patients having different values of lung compliance and pressure losses in the trachea;  to study different architectures by using different components or by changing their position in the circuit.