Lung ventilation in Simulink

Pulmonary ventilation is a critical operation both during surgery and for assisted or replacement ventilation of the patient.

This report proposes the analysis of a lung mechanics model through an electrical analogy with an equivalent circuit consisting of resistors and capacitors. This electrical-mechanical analogy makes it possible to analyze how the variation of lung parameters affects ventilation in the case of pressure-controlled ventilators and therefore to simulate the response of the pulmonary system to ventilation.

A mathematical model is also proposed for the implementation of a pressure control ventilator with square wave, within the Simulink modeling environment, with a graphical interface and user-selectable parameters.

Various pathological conditions are then simulated by going to see how the variation in lung compliance has the predominant effect. Small variations in the capacity values ​​are reflected in large variations in the volume of incoming air up to dangerous states for the patient himself.

The following report is an extract from a project carried out for the course of Modeling and Simulation of Physiological Systems, held for Medical Engineering at the University of Rome Tor Vergata.

Introduction

The evaluation of the respiratory system with electrical analogies makes it possible to analyze different distinctive clinical pictures of pulmonary pathophysiology.

There are many models in the literature and several authors attribute to the pulmonary system different inertial and elastic contributions often translatable into an electrical model of an RLC circuit [1]. In this report the inertial contributions are neglected and an RC model is analyzed which allows to describe the pulmonary system with good approximation.

The introduction of these models in simulation software allows you to analyze and model different respiratory behaviors. A model of respiratory mechanics can be integrated with the modeling of a pressure-controlled ventilation system by imposing the waveform of pressure entering the airways. The combination of the two simulation models makes it possible to study the diagnosis, evolution and treatment of particular pathologies and to analyze how the responses of the respiratory system vary when the parameters of the respirator vary.

The need for ventilation support is a frequent cause of patient transfer to the ICU and analyzing complex lung physiology with simplified models would lead to better management of the patient and ventilation devices.

Furthermore, a digital simulation system allows to better analyze the flow curves towards which there is more and more interest in the clinic [2]. These curves contain a lot of information on respiratory mechanics, patient effort and ventilation mode, including parameters. Correctly analyzing these curves allows us to understand if the patient’s breathing pattern is synchronized with the machine or if too much effort is required from the patient which could require prolonged ventilation or lead to serious lung damage.

Background

Correctly simulating the patient requires the choice of a sufficiently accurate respiratory mechanics model and then the fitting of this model to the patient’s peculiarities.

The following section analyzes a model for respiratory mechanics and compares some parameters of resistance and compliance highlighted in the literature as average properties of the population.

This model is implemented in Simulink together with different pressure waveforms. Subsequently, the description of a pulmonary ventilator model in Simulink is also dealt with in order to allow the user to customize the ventilation parameters.

Circuit analogy

The pulmonary circuit can be analyzed by making an analogy with the electrical circuits.

In particular, it is possible to make a parallelism between the flow of air and the electric current (flow of charges) and the pressure as the presence of an electric potential. The mechanical resistance is then reviewed as the ratio between the pressure increase with respect to the flow, analogous to the electrical resistance. Similarly, compliance is nothing more than the relationship between the increase in volume and the increase in pressure, in electrical analogy it is a condenser.

The system in fig. 1a is a model of respiratory mechanics that neglects the presence of inertial contributions (there are no inductances) and considers the presence of two compartments. The upper airways are separated, with their resistive contribution 𝑅𝐶 from the lower airways 𝑅𝑃 . The two compartments are in series with each other and in series with the air tanks, i.e. the capacities representing the compliance contribution of the 𝐶𝑊 wall and of the lung 𝐶𝐿 . These contributions are in series precisely because the throughput is the same.

Added to this is the shunt capacity 𝐶𝑆 which takes into account various contributions such as anatomical dead space, airway deformability and air compressibility. This capacity is parallel to the pulmonary path, representing the air flow that does not manage to carry out gas exchanges. Normally this volume is very small in physiological respiratory conditions and at low respiratory frequencies, it is negligible.

The pressures in the nodes are also identified. Airway pressure 𝑃𝑎𝑤 , pleural pressure 𝑃𝑝𝑙 and alveolar pressure 𝑃𝐴 . Clearly, the entrance to the system, given by the mouth and nasal cavities, is represented by the pressure at the opening of the airways 𝑃𝑎𝑂 .

System response

The circuit in fig. 1a can be described by the following equations:

\begin{equation}
\left\{\begin{array}{l}
P_{a O}=Q R_{C}+\frac{1}{C_{S}} \int\left(Q-Q_{A}\right) \\
\frac{1}{C_{s}} \int\left(Q-Q_{A}\right)=Q_{A} R_{P}+\left(\frac{1}{C_{L}}+\frac{1}{C_{W}}\right) \int Q_{A}
\end{array}\right.
\end{equation}

We then obtain the system transfer function:

\begin{equation}
\begin{aligned}
H(s) &=\frac{Q(s)}{P_{a O}(s)}=\\
&=\frac{s^{2}+s \frac{1}{R_{P}}\left(\frac{1}{C_{S}}+\frac{1}{C_{e q}}\right)}{s^{2}\left(R_{C}\right)+s\left(\frac{R_{C}+R_{P}+\frac{R_{C} C_{S}}{C_{e q}}}{C_{S} R_{P}}\right)+\frac{1}{C_{e q} C_{S} R_{P}}}
\end{aligned}
\end{equation}

Where series of capacity is expressed as:

\begin{equation}
\frac{1}{C_{e q}}=\frac{1}{C_{L}}+\frac{1}{C_{W}}
\end{equation}

System properties

Clearly, the solution of this problem requires the knowledge of the patient’s own lung mechanics so as to be able to fix the different coefficients of the model.

It is therefore assumed that the patient has normal lung mechanics and the numerical coefficients are selected by Khoo [3] and are reported in Tab. 1.

System modeling

Solving the system requires solving an ODE and the approach, when the equations become complex, is to transfer the model to a computer. The classic approach is numerical, through the use of a numerical code for solving differential equations.

However, there is the possibility of using Simulink.

Simulink

Simulink [4] is software developed by MathWorks that provides a graphical approach based on a modeling environment that allows the user to convert the problem into a network of mathematical function blocks. Furthermore it could be integrated with the Matlab language and related programming functions.

A first synthetic approach could be to diagram an input-output system through the transfer function in eq. (2). However, this model would be too synthetic and would not allow access to some internal variables, such as individual flows.

We then choose to model the complete system in eq. (1). This system is modeled in the lungs block. A subsystem is then added to simulate a pulmonary ventilator, i.e. the input as 𝑃𝑎𝑂 , and a block to view and save the data. A general high-level diagram is present in the appendix (fig. 19).

Ventilator subsystem

In the fan subsystem, the goal is to provide a 𝑃𝑎𝑂 with a precise waveform. Different waveforms are considered and the frequency itself, in breaths per minute, can be varied.

To do this Matlab function block can be used to define a Matlab code containing the instructions to generate the pressure waveform as the input time varies. More information can be found in the appendix.

This block also requires two auxiliary variables through which you can choose the respiratory rate and the type of waveform directly from the graphic interface of the block mask (fig. 2b). A callback function is also set through the block mask that allows you to update the name of the output file, as soon as the parameters are changed, with a structure like:

"wave \:form + frequency + .mat"

Lung subsystems

The descriptive equations of the system (eq. (1)) and the related circuit segments can be represented directly in a Simulink model, presented in fig. 3.

By modeling in Simulink it is possible to add the signal contributions (representing segments of the circuit or, equivalently, members of the equation), multiply by a constant by applying a gain to the signal, derive and integrate. Clearly, to go from flow to volume it is sufficient to integrate over time.

There are also 3 scope-type blocks to display the waveforms directly within the simulation, the input block (takes the 𝑃 signal and the 3 output blocks are used to save the data.

Results

Ideal wave form for ventilation

A simple approach to analyzing lung mechanics is to ventilate with an ideal waveform. A sinusoidal shape with an amplitude of 2.5 cmH2O (peak-to-peak amplitude of 5 cmH2O) with a frequency of 15 breaths per minute [3], similar to resting breathing, is used as the first analysis.

From the graphs in fig. 4 it is possible to see how the volume follows a similar trend, albeit slightly out of phase. This indicates that, at low frequency, the trend is dominated by the contribution of compliance (where there is proportionality with the integral of the flow). The peak volume amplitude reaches 0.5 L and the flow around 0.7 L / s. The flow trend, on the other hand, is significantly out of phase with respect to the input pressure contribution.

Frequency dependence

By increasing the frequency, up to triple the number of breaths per minute (fig. 5), the lung appears more rigid and in fact the amplitude of the volume tends to decrease despite the increase in flow. The peak amplitude of the volume goes from 0.5 L to 0.4 L although the flow has gone from 0.7 L / s to almost 2 L / s. This represents how although the air exchange is greater, the lung is forced to expand less in order to increase the respiratory rate, limiting the overall intake of air.

Furthermore, the volume trend tends to be more out of phase than the pressure trend while the flow seems to be more in phase, a symbol of a greater resistive contribution (direct pressure-flow proportionality).

Therefore, while at low frequency (rest) the contribution is more related to compliance, at high frequency the resistive behavior seems to be more present.

Square wave

An alternative waveform is the square wave. A waveform with a peak amplitude equivalent to the previous one can be set but a different system response is observed. In fig. 6 we see how the flow quickly returns to zero and a volume plateau phase begins. Then begins the inspiratory phase where the flow presents a symmetrical but negative curve with a peak of the same amplitude and the volume gradually increases.

Furthermore, by varying the frequency the behavior of the flow and volume curves does not change but the plateau phase of the volume is reduced. Clearly, these waveforms are analytical models that in real application would present various problems. In clinical practice, negative pressures are never used but a positive pressure gradient is always maintained, so these curves are to be excluded for realistic application cases.

Also note how the very presence of a negative volume, in the sinusoidal waveform, is not desirable for a ventilator whose purpose is to help the patient breathe by reducing the lung’s respiratory work or replacing it completely.

The previous analyzes are useful for understanding the mathematical model and for having a basic idea of respiratory mechanics. However, in the clinic, the approach to ventilation is slightly different.

Pressure-controlled ventilation

In order to create a digital simulator of a pulmonary ventilator, a waveform typically diffused in the clinic for pressure-controlled ventilation is considered [5].

Assisted pulmonary ventilation is divided into two macro categories: volume control and pressure control. There are several advantages / disadvantages to both and the choice is based on several specific patient considerations. In pressure-controlled ventilation (VCP) you set a specific target, the airway pressure 𝑃𝑎𝑂 and the ventilator must maintain it. Depending on the machinery used, some parameters may be varied such as the percentage of oxygen, the tidal volume or ventilation per minute, the respiratory rate, the inspiratory time, the flow or set limit values ​​for the pressure [6].

A generic curve, typical of VCP, consists of a square wave with a positive offset which is called PEEP, positive end-expiratory pressure. The period 𝑇 is composed of a breathing time, inspiration time and an inspiration time. Optionally, a certain rise time may be present, expressed as a function of the inspiration time [5].

This can be expressed as a piecewise linear curve (fig.7b) of period 𝑇 = 𝑇𝑖𝑛𝑠𝑝 + 𝑇𝑒𝑠𝑝 , mathematically described as :

\begin{equation}
P(t)=\left\{\begin{array}{cl}
P_{a O} \cdot \frac{t}{\tau}+\text { PEEP } & 0 \leq t<\tau \\
P_{a O}+\text{ PEEP } & \tau \leq t \leq T_{in s p} \\
\text { PEEP } & T_{i n s p} \leq t
\end{array}\right.
\end{equation}

PEEP generally remains between 5 ÷ 10 cmH2O keeping the volume in the range of 4 ÷ 6 mL / kg of body weight. The 𝑃𝑎𝑂 is kept below 35 cm H2O [7, 8].

GUI

In this regard, the graphic interface (fig.2b) is also modified in such a way that it allows the selection of a square wave waveform with a rise time and a pressure offset, leaving the user the possibility to choose frequency, inspiration time (as a function of the period), rise time (as a function of inspiration time), 𝑃𝑎𝑂 and PEEP. The complete interface is present in the appendix (fig. 18).

Once set, this simulator allows you to calibrate the ventilation and to train in the use of a real ventilator by going to see how the different ventilation parameters affect the respiratory mechanics and above all how to modify the ventilation according to the patient’s response.

Physiological breathing

In the following sections the ventilation parameters [9] are set at a rate of 18 breaths / min, an inspiration time of 30% and a rise time of 15%, with a PEEP of 5 cm H2O and a 𝑃𝑎𝑂 = 25 cm H2O These parameters allow to approximate the waveforms typically used in the clinic [5].

In order to satisfy the lung response, in terms of tidal volume in order not to bring the patient into dangerous conditions [10], it is necessary to satisfy an overall lung compliance of approximately:

\begin{equation}
C_{T}=\left(\frac{1}{C_{W}}+\frac{1}{C_{L}}\right)^{-1} \approx 10 \mathrm{~mL} / \mathrm{cmH} 2 \mathrm{O}
\end{equation}

Then follows the analysis considering the average parameters of a lower order of magnitude: 𝐶𝑤 = 0.02 L / cmH2O and 𝐶𝐿 = 0.02 L / cmH2O. It is then analyzed how the ventilation varies according to the variation of the patient’s parameters.

Airway resistance

There are several pathologies and clinical pictures that result in an increase in resistance of the internal airways. An increase in resistance can be caused by histological alterations, in the alveolar geometry or in the alteration of the air-liquid interface.

The muscular work that must be generated by the respiratory muscles is dependent on the elastic and resistive properties of the respiratory system for the entire range covered by the volume variation. Normally the lung has a compliance of about 200 mL / cmH2O, i.e. when the transpulmonary pressure (difference between the pressure in the alveoli and the pleural one) increases by 1 cmH2O then the lung expands by 0.2 L. Clearly, the smaller the compliance more the slope in the pressure-volume diagram will be large.

Pathophysiology

The elastic force is the greatest contribution (> 2/3 of the total) in causing the lung to collapse due to the surface tension inside the alveoli. Although this is not a problem in the healthy lung, it becomes a problem in cases where the amount of pulmonary surfractant is reduced, such as in acute respiratory distress syndrome (ARDS) [11], where compliance decreases.

Emphysema is characterized by an increase in respiratory resistance and a decrease in lung compliance [11].

There are also several factors not related to specific pathologies that tend to increase resistance such as the presence of obstructions in the endotracheal tube, cough, secretions, bronchospasms and accelerated respiratory rhythm [6].

Other influencing factors

Added to this are external factors that can vary both resistance and compliance such as restrictions on the ventilation system, abdominal contraction or increased abdominal pressure or damage or deformation of the chest wall.

It can be seen how, by varying both resistances present, the volume response is almost zero. Instead, the flow varies slightly from the peak reached. From figs. 9 to 12 we see how the flow is more sensitive to the variation of 𝑅𝑐 .

Remembering that the tests are performed at the same pressure, an increase in resistance will result in a reduction in flow.

Pulmonary compliance

Lung compliance is also affected by some pathologies.

In chronic obstructive pulmonary disease (COPD) there is a reduction of the pulmonary compliance. In the case of cardiogenic (CPE) or non-cardiogenic (ARDS) pulmonary edema there is a reduction in compliance even up to 0.044 L / cmH2 O and 0.035 L / cmH2 O, respectively [11]. To this are added other factors that contribute to reducing lung compliance such as atelectasis, pneumothorax, displacement of the endotracheal tube and pneumonia [6].

The variation in compliance greatly affects the volume of air entering the lungs. In particular, since 𝐶𝑠 is typically very small and low, it is possible to neglect this contribution and therefore the quantity of flow that is directed on the shunt portion and does not enter the lungs is very low, it is possible to neglect this contribution and consider the volume entering the lungs as the volume delivered by the ventilator.

Shunt capacity

However, this is not true if 𝐶𝑠 increases. In fact, as can be seen from fig. 16, the more 𝐶𝑠 increases, the more the flow on the shunt branch increases and the greater the gap between the total flow (which instead remains almost constant with an average deviation of less than 10−3 [L / s]) and the alveolar flow.

In the figure fig. 16b shows the norm of the difference between the flow delivered by the respirator and the alveolar flow, normalized with respect to the maximum value.

\begin{equation}
\Delta \mathrm{Flux}=\frac{\| \text { flusso su } C_{s}-\text { flusso alveolare } \|}{\max \left(\text { flusso su } C_{s}-\text { flusso alveolare }\right)}
\end{equation}

Alveolar flow

To be sure of the amount of air actually entering the lungs and capable of gas exchange, it is necessary to read the flow value on the lung branch, or rather on the capacity 𝐶𝑇. Comparing the alveolar flow values fig. 17, as the system parameters vary, with the respective variations of the total flow it is possible to identify how much flow actually enters the lung.

In the case of the variation of the resistances, the difference between the flows remains very low compared to the variation of the 𝐶𝑇 for which the variations of the flow are seen to correspond almost completely. On the other hand, the shunt capacity is different for which there is a direct increase in the shunt flow to the detriment of the alveolar flow (fig. 16).

Conclusion

The results show that it is essential to correctly set both the ventilation parameters and to estimate the patient’s behavior.

Variations in the circuit parameters can alter the behavior and quantity of air required from the ventilator. In particular, the term with greater importance is given by the series of lung capacities whose increase leads to an excessive increase in the input volumes up to values ​​that can be dangerous for the patient.

The increase in shunt capacity tends to increase the diverted flow and therefore to remove air from the lung, despite the fact that the fan still introduces more air. Furthermore, the variation of the resistances, in the ranges considered, is not reflected so much on the variation of the overall flow as on the slope of its waveform.

So, the availability of a pulmonary ventilator model coupled to a respiratory mechanics model allows you to simulate the ventilator-pathology interaction and analyze the best method of operating depending on the parameter that is deviated from the pathology.

Code availability

All the code is available in the GitHub repository.

References

• [1] Pardis Ghafarian, Hamidreza Jamaati, and Seyed Mohammadreza Hashemian. “A Re- view on Human Respiratory Modeling”.
• [2]  Natsumi Hamahata, Ryota Sato, and Ehab Daoud. “Go with the flow—clinical im- portance of flow curves during mechani- cal ventilation: A narrative review”.
• [3]  Michael C. K. Khoo. Physiological control systems: analysis, simulation, and estima- tion. en. Second editon. IEEE Press series in biomedical engineering.
• [4]  MathWorks. “Simulink”. In: (2022). URL:https://it.mathworks.com/products/ simulink.html. Noman Q. Al-Naggar. “Modelling and Sim- ulation of Pressure Controlled Mechanical Ventilation System”. en. In: Journal of Biomedical Science and Engineering 08.10
• [6]  Irene Grossbach, Linda Chlan, and Mary Fran Tracy. “Overview of Mechanical Ven- tilatory Support and Management of Patient- and Ventilator-Related Responses”. en. In: Critical Care Nurse 
• [7]  Lorenzo Ball, Maddalena Dameri, and Paolo Pelosi. “Modes of mechanical ventilation for the operating room”. en. In: Best Practice & Research Clinical Anaesthesiology 29.3
• [8]  Dean Hess. “Ventilator waveforms and the physiology of pressure support ventilation”. In: Respiratory care 50
• [9]  R. Duncan Hite. “Modes of Mechanical Ven- tilation”. en. In: A Practical Guide to Me- chanical Ventilation. Ed. by Jonathon D. Truwit and Scott K. Epstein.
• [10]  (Roy G. Browe. “Ventilation with Lower Tidal Volumes as Compared with Tradi- tional Tidal Volumes for Acute Lung In- jury and the Acute Respiratory Distress Syn- drome”. en. In: New England Journal of Medicine 342.18
• [11] J. Milic-Emili et al., eds. Basics of Respira- tory Mechanics and Artificial Ventilation. en. Milano: Springer Milan, 1999. 

Appendix

Code for generating the inlet pressure

function y = fcn(t,flag_ventilation,breath_for_minute,tau_int,Paw,PEEP,T_insp) % % % to  avoid error use the ‘odeN‘ solver !
 f=breath_for_minute/60; % breathing frequence in Hz
 % pressure diffrence from maximus value and PEEP
 DeltaPressure=Paw-PEEP;
 A=0; % initialize output variable
 if flag_ventilation==1 % flag=1 for sine wave
  % sine wave from Khoo Physiological Control Sistem % amplitude peak -peak of 5 cmH2O
  amplitude =5;
  A=amplitude*sin(2*pi*f*t);
 end
 if flag_ventilation==2 % flag=2 for perfect square wave
   dt=rem(t,(1/f)); % current time (of single period)
   period=1/f;
   half_period=period/2; % duty cycle of 50 % by selecting half period if   dt<half_period % semi-period where the output is on
   amplitude=5; % amplitude of 5 [cmH2O]
   A=amplitude; % set the output constant end
   if dt>= half_period % semi-period where the output is zero 
    A=0; % set the output to zero
   end 
  end
 if flag_ventilation==3 % flag=3 for square wave with rise time tau
  period=1/f;
  dt=rem(t,(1/f)); % current time (of single period)
% the inspiratory time (where the output is on) is defined usign percentage of total period with T_insp from ventilator GUI % also the rise time is defined using percentage of inspiratory time with tau_int from ventilator   GUI insp_Time=(T_insp/100)*period;
  tau=(tau_int/100)*insp_Time; if dt < tau % rise time
% linear growth of the pressur with offset
  A=DeltaPressure*(dt/tau)+PEEP;
  end
  if (dt>=tau) && (dt<insp_Time) % the pressure is max (square wave on)
  A=DeltaPressure+PEEP;
  end
  if dt>= insp_Time %the pressure is min (square wave off)
  A=PEEP;
  end 
 end
 y=A; % set the pressure wave as function output 
end

Complete graphic interface

High-level model