Towards a multi-scale computer modeling workflow for simulation of pulmonary ventilation in advanced COVID-19.

Physics-based multi-scale in silico models offer an excellent opportunity to study the effects of heterogeneous tissue damage on airflow and pressure distributions in COVID-19-afflicted lungs. The main objective of this study is to develop a computational modeling workflow, coupling airflow and tissue mechanics as the first step towards a virtual hypothesis-testing platform for studying injury mechanics of COVID-19-afflicted lungs. We developed a CT-based modeling approach to simulate the regional changes in lung dynamics associated with heterogeneous subject-specific COVID-19-induced damage patterns in the parenchyma. Furthermore, we investigated the effect of various levels of inflammation in a meso-scale acinar mechanics model on global lung dynamics. Our simulation results showed that as the severity of damage in the patient's right lower, left lower, and to some extent in the right upper lobe increased, ventilation was redistributed to the least injured right middle and left upper lobes. Furthermore, our multi-scale model reasonably simulated a decrease in overall tidal volume as the level of tissue injury and surfactant loss in the meso-scale acinar mechanics model was increased. This study presents a major step towards multi-scale computational modeling workflows capable of simulating the effect of subject-specific heterogenous COVID-19-induced lung damage on ventilation dynamics.

Introduction

Coronavirus Disease 2019 (COVID-19) infection due to the Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) virus can cause extensive damage to many tissues and organs, including the lungs. Millions of lives have been lost to the disease, while others face long-term effects of this viral infection [1,2]. Though computer models may assist with managing and tracking the spread of COVID-19, yet new variants of the SARS-CoV-2 virus may spread more easily and can pose an increased risk of severe disease and other complications linked to COVID-19 [3,4]. Alveolar damage is seen in imaging and histopathological studies of COVID-19 infected lungs [5], and heterogeneous damage throughout the acinar regions of the lung is observed in computed tomography (CT) images [5]. In particular, CT imaging can provide useful spatial information on patterns of lung injury, including ground-glass opacities (GGO) and areas of consolidation [6]. Also, some patients suffering from COVID-19 acute respiratory distress syndrome (CARDS) exhibit increased hypoxemia compared to typical acute respiratory distress syndrome (ARDS) [7]. However, there is currently little quantitative information on how acinar level patho-mechanics contribute to whole lung function in COVID-19 patients. A better understanding of how COVID-19 impairs regional ventilation may give insight into overall lung dynamics specific to CARDS.

To provide a four-dimensional view of airflow patterns in the lung, CT images can be used to develop geometric models of the lung, which then can be combined with fluid flow and tissue mechanics physics-based models. Such physics-based computer models of the lung present an opportunity for gaining valuable insights into pulmonary ventilation dynamics [8,9]. For instance, computational modeling of lung dynamics across multiple spatial scales may allow a deeper understanding of how mechanical changes at the alveolar and acinar levels affect lobar and whole-lung dynamics in CARDS. Previous physics-based computer models of the lung have provided a detailed four-dimensional view of pulmonary ventilation in both healthy and disease states. In the aforementioned computer models, CT images were used to determine lobar volumes and airway branching patterns of large airways with small peripheral airways constructed utilizing volume-filling tree-generating algorithms [10,11] and the coupling of airway trees to compliant acinar regions provided realistic flow distribution among the lung lobes [12]. Along with realistic geometries for airways and compliant acinar regions, in silico models have also incorporated other factors such as tissue deformation, gravity, acinar-level interdependence, and surfactant all of which contribute to distribution of ventilation in the healthy lung [[13], [14], [15]].

Disease can lead to remodeling of airways and parenchymal tissue thus inducing changes in airflow patterns and tissue mechanics that can be implemented in in silico models [16]. Insight into the effects of disease on ventilation and pressure distribution of the lung has been gained through modification of model geometry and mechanics in previous studies [17,18]. Additionally, information on distribution and manifestation of damage throughout the lung is important for realistic representation of disease states [18]. Registration of high resolution and 4D CT images has been used to identify damaged areas of parenchyma and changes in airway geometry and regional ventilation in disease states [16,19]. Utilization of imaging techniques also opened the door for patient-specific modeling of airflow and pressure distribution [10,18]. Multi-scale in silico lung models with geometries resolved from CT imaging have proven useful in understanding various pulmonary diseases like cystic fibrosis, chronic obstructive pulmonary disease (COPD) and asthma [18,20,21]. However, there is a need for computer modeling studies on pulmonary ventilation dynamics in COVID-19 patients..

The main objective of this study is to develop a physics-based in silico modeling workflow for studying the pulmonary ventilation of COVID-19-infected lungs by bringing existing methodologies for coupling airflow and lung tissue mechanics together [10,11,22]. The presented in silico modeling approach is the foundation and first step towards a virtual hypothesis testing platform for a better understanding of COVID-19 pulmonary dynamics utilizing 4D CT data from COVID-19-infected lungs and accoutning for patient-specific lung geometry and disease distribution with varying levels of damage. In addition, we aimed to develop an in silico multi-scale approach to simulate and compare the regional changes in lung dynamics associated with heterogeneous subject-specific COVID-19-induced damage in the parenchyma of the lung. To this end, we present an investigation of the effect of various levels of inflammation in a meso-scale acinar mechanics model on global lung dynamics.

Methods

Imaging

This study utilized the 4D CT scan of a male patient recently hospitalized in Vidant Medical Center (Greenville, North Carolina, USA) for an advanced case of COVID-19. The 4D CT scan was obtained during tidal breathing using an Optima CT580 RT scanner (GE Healthcare, Waukesha, WI). The methodology used in this paper was approved by the East Carolina University and Medical Center Institutional Review Board (UMCIRB) with study ID 20–001447. Informed consent was obtained from the patient. Sorting of the CT images into phases of the breathing cycle was accomplished using the Varian 4DCT Real-time Position Management (RPM) system (Varian Medical Systems, Palo Alto, CA). When a series of CT images was attained over a time comparable to that of a normal breathing cycle, the Varian RPM camera captured the patient's real-time external chest motion amplitude, and Advantage 4D (GE Healthcare, Waukesha, WI) software retrospectively sorted the CT data into corresponding phases of the respiratory cycle from 0% to 90%, with 0% corresponding to end-inspiration and 50% corresponding to end-expiration [23]. The images were taken at 120 kVp, 300 mAs, and 20 mm collimation and were reconstructed through a 512 × 512 matrix with a 2.5 mm slice thickness and reconstructed retrospectively to a slice thickness of 1.25 mm.

Segmentation

Following imaging, the geometry of the major airways visible in CT and lungs lobes was segmented. The major conducting airway geometry was segmented from the end-inspiratory phase using a combination of dynamic region growing and manual editing in Materialise Mimics 23.0 (Materialise NV, Belgium) (Fig. 1 a). Centerline detection and extraction as described in Bordas et al.‘s work [10] for the major airways were also achieved using Mimics 23.0. Five different lung lobe (right upper, right middle, right lower, left upper, and left lower) geometries were segmented at the end-inspiratory phase and end-expiratory phase using the Chest Imaging Platform [24] available in 3D Slicer [25,26]. Total segmented lobe volume was validated against total segmented lung volume. 3D Slicer was used to segment the COVID-19 affected GGO and consolidated regions of the lungs at end-inspiratory and end-expiratory phases using the LungCTAnalyzer [27] extension of the Chest Imaging Platform (Fig. 1b). This segmentation was based on Hounsfield unit (HU) values in the CT images. Inflated lung Hounsfield unit thresholds were determined based on the study by Kassin et al. [28] with a range of −1000 to −650 HU. The difference between ground glass opacities and consolidated regions were determined based on the 3D Slicer Chest Imaging Platform [24] and Lung CT Analyzer [27] (https://github.com/rbumm/SlicerLungCTAnalyzer) extensions preset values for COVID-19 lung analysis.

Fig. 1

CT image as segmented in 3D Slicer and Mimics software. (a) Purple regions show major airways at end-inspiration. (b) Areas affected by COVID-19 damage, such as GGO (orange) and consolidated (blue) are highlighted based on HU values used for segmentation. The remaining dark grey areas are aerated lung tissue; (c) Geometrical representation of the entire generated tree, viewed from the front, including the airways generated using the space-filling algorithm. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Geometry

Segmented major airway geometry was limited to the first four to six generations. Further airway segmentation was constrained by CT image resolution. Subsequently, a space-filling airway generation algorithm with random heterogeneity based on the work of Tawhai et al. [11,29] was used to create up to 16 generations of conducting airways. Airway diameter, length, and branching angles were based on the segmented geometry and lung lobe models (Fig. 1c) following the methods described in Refs. [[29], [30], [31]]. The number of acini per lung was approximately 15,000 [32]. The lumen diameter for each one-dimensional line segment was assigned based on the Horsfield number [10,11]:where , , , represented the current Horsfield order, the airway diameter, the maximum Horsfield order and the maximum diameter, respectively. represented the anti-log slope of airway diameter plotted against Horsfield order and was assigned to be 1.15.

Airflow and acinar mechanics model

Here, our multi-scale computational models of the lung are developed through the C++ simulation package CHASTE (Cardiac, Heart, and Soft Tissue Environment) [10,33]. Airflow in the airway tree geometry was described by a reduced dimensional airway model implemented in CHASTE, which was coupled to the tissue mechanics acinar models [10,34]. The flow was presented as a modified Poiseuille flow following the approach developed by Swan et al. [13] and Ismail et al. [12] and by assuming isotropic expansion of acini. Corrections to the dynamic resistance based on work by Pedley et al. [10,35] were applied as shown in Equation (2):where R is the dynamic resistance, R is the Poiseuille resistance, D and l are the diameter and length of an airway, respectively, Re is the Reynold's number, and was set to be generation-dependent based on the work of van Ertbruggen et al. and Ismail et al. [12,36].

Atmospheric pressure was assigned at the trachea, hence airflow into the lungs was driven by variations in the volume of acini as a function of changes in the transpulmonary pressure during breathing. All nodal pressures and edge fluxes were solved for simultaneously using multifrontal lower–upper factorization solver UMFPACK. Flow from the airway model was fed into each acinar model to calculate the change in acinar volume using the stretch ratio during each time step taken by the solver.

A sigmoidal acinar mechanics model based on the work of Fujioka et al. and Venegas et al. [37,38] was coupled to the generated airway tree in the simulations (Fig. 2 ). In this model, as shown in Equation (3), V is the acinar volume, is the transpulmonary pressure for each acinus (defined as the difference between pressure in the acinus and the pleural pressure), and A, B, C, and D are constants that vary based on surfactant level and consequently tissue compliance:

Fig. 2

The sigmoidal model derived from Refs. [22,39] used in a simulation for healthy lung and COVID-19 affected lungs. The pressure – volume curve for a healthy lung is shown in blue, with pressure being equal to the transpulmonary pressure. As damage progresses in the lungs, there is a progressive reduction of surfactant amount and compliance of the acinar units. The decrease in surfactant shifts the pressure-volume curves to the right, as seen in the 20% reduced surfactant (orange), 40% reduced surfactant (grey), and 60% reduced surfactant (yellow) cases. The dashed grey lines indicate the pressure range used for the simulation. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Time-derivative of equation (3) was solved at the end of each terminal bronchiole, where an acinus was connected to an airway, thus coupling the pressure in the airway tree and the acinar model.

Boundary conditions and simulation settings

Boundary and loading conditions

Pulmonary ventilation dynamics during tidal breathing was simulated in accordance with the aforementioned coupled airway-acinar model. To simulate tidal breathing, a varying pleural pressure was applied at the acini while a constant atmospheric pressure boundary condition was applied at the trachea. The varying pleural pressure was assumed to be [12,40]:where (-490 Pa), (-314 Pa), and the phase shift .

Simulation of COVID-19-afflicted lungs

The coupled airway-acinar model was used to simulate COVID-19 lung damage based on the segmented CT images and the associated region-specific levels of damage, corresponding to GGO and consolidated regions. Since the CT images were obtained from a patient recovering from advanced COVID-19, hypothetical healthy lung simulations were also performed to provide a basis for comparison of the results. To simulate the changes between healthy and COVID-19 afflicted lung function, different mechanical behavior of the sigmoidal acinar model (Fig. 2) was implemented based on the amount of surfactant and compliance of the acinar units. Simulations were run for the healthy case considering the normal amount of surfactant. The COVID-19 affected lung simulations considered the reduced amount of surfactant and decreased compliance in both GGOs and consolidation regions [[41], [42], [43]].

Previously segmented GGO and consolidated regions in lung CT images were used to identify acinar regions affected by COVID-19. We assumed GGOs to represent partial filling of air spaces while consolidated regions may signify more severe damage where a large proportion of airspace is filled with liquid and inflammatory infiltrates [5,6]. While the airspace becomes compromised and inflamed, as seen in CT images of GGO and consolidation regions, the amount of surfactant and compliance of the acini can be altered [38,44,45]. Different acinar properties were applied to the GGO and consolidated regions to simulate these varying levels of tissue damage. Two simulations were run to visualize the effects of varying levels of COVID-19 severity: for the first simulation, a 20% reduction in the surfactant amount was used for the GGO regions, and 40% reduction in the surfactant for the consolidation regions was considered. For the second simulation, a 20% reduction of surfactant in the GGO regions and a 60% reduction in the consolidated regions was considered (Fig. 2).

Coefficients for the healthy sigmoidal model were fit to physiologically relevant coefficients based on patient-specific lung volumes [39]. Lung volumes were estimated based on the gender, height, and age of the COVID-19 positive patient used for CT scans in this study: male, 167.6 cm, 51 years, and 88.5 kg, respectively [46].

The following equations from Boren et al. [47] were used to estimate lung capacities and volumes for a healthy adult male: where TLC is the total lung capacity in liters, H is the subject's height in centimeters, FRC is functional residual capacity in liters, RV is residual volume in liters, Ag is the subject's age in years, and VC is vital capacity in liters. Based on these equations for the patient in this study, we found that TLC = 5.77 L, FRC = 2.42 L, RV = 1.53 L, and VC = 3.99 L. Estimated lung volumes were then used to define coefficients A and B, as A is approximated by the residual volume and B is approximated by the vital capacity [39]. The values of C and D in the sigmoidal model correspond to the inflection point of the curve and the pressure range in which the volume change takes place, respectively [39]. As such, C and D are not as easily determined by available physiological data and were estimated such that in the healthy case, they were estimated based on the calculated FRC and the tidal volume pressure range.

Then, to define the regions of GGO and consolidation, C and D coefficients signifying reduced surfactant as reported by Fujioka et al. [38], were scaled to fit our patient data. Fujioka et al.‘s [38] coefficients were limited to normal surfactant, 20% reduced, and 40% reduced surfactant amounts. Coefficients for 60% reduced surfactant was not directly available. Hence, the coefficients for 60% reduced surfactant were extrapolated. Assumptions for extrapolation were based on work by Fujioka et al. [38] that only coefficient C varied markedly when surfactant level was changed (Table 1 ).

Table 1

Coefficients for use in the sigmoidal model of equation (4), defining the acinar unit's pressure-volume relationship.

	A (L)	B (L)	C (Pa)	D (Pa)
Healthy	1.53	4.24	1078.73	449.14
20% reduced surfactant	1.53	4.24	1420.00	451.11
40% reduced surfactant	1.53	4.24	1818.15	356.96
60% reduced surfactant	1.53	4.24	2359.48	356.96

Three different simulations were executed for this study:

A hypothetical control study where the lung was assumed to be healthy; GGO and consolidated regions are considered to be normal inflated regions with normal surfactant levels.

COVID-19-afflicted lung with 20% reduced surfactant for the GGO region and 40% reduced surfactant for the consolidated region

COVID-19-afflicted lung with 20% reduced surfactant for the GGO region and 60% reduced surfactant for the consolidated region.

The total simulation time was 12 s, with each breathing cycle assumed to be 4 s between consecutive end-inspirations. The time step used by the solver was 0.001 s, and data were saved every 100 time steps. A smaller 0.0001 s time step was also tested for the healthy case and showed no significant difference in results other than an increase in simulation time. Three breathing cycles were executed, with the first two cycles discarded so that only steady-state conditions were analyzed. Following the onset of steady-state conditions, the data from that breathing cycle was used to calculate the flow rate at the trachea, flow rate into the individual lobes, total tidal volume of the lungs at different time points, the flow rate through individual lobes. Flowcharts summarizing the entire model development and simulation process are shown in Fig. 3, Fig. 4 , respectively.

Fig. 3

Flowchart for generating the airway tree geometric model from CT images. The CT image was segmented to determine major conducting airways, lungs, and lobes as well as COVID-19-affected regions. CHASTE [33] was then used to generate the complete airway tree down to terminal bronchioles.

Fig. 4

Flowchart for running the tidal breathing simulation for healthy and diseased lungs. The airway tree model was coupled with the sigmoidal acinar model in CHASTE [33] to simulate tidal breathing. Different levels of surfactant reduction were applied to the acinar model to simulate lung function in disease states.

Results

In total, three simulations were performed consisting of a hypothetical healthy lung with normal surfactant levels as well as two COVID-19-afflicted simulations created through reduction of surfactant levels based on the degree of damage present in the CT images. The first simulation will be referred to as the hypothetical “healthy” case, while the second and third simulations will be referred to as “diseased 20–40” and “diseased 20–60” cases, respectively.

Airflow during the simulations was generated due to the negative pressure as a result of the acinar expansion. Fig. 5 shows the flow rate of air at the trachea and into each lobe after reaching steady-state over an entire breathing cycle. Different colored lines represent each simulation scenario, with the highest flow rate in the healthy simulation and the lowest in the diseased 20–60 simulation. The lobar flow rate values were determined by calculating the flow rate at the first airway branch that enters each lobe.

Fig. 5

Flow through the trachea (a) and each lung lobe (b–f) during tidal breathing over one breathing cycle, plotted by tracking the flow every 100 time steps.

The flow rate at various points in the lung can be integrated to determine the total and lobar tidal volumes (Fig. 6 ). The maximum values of tidal volumes for the whole lung and each lobe are also reported in Table 2 , in addition to the percent change in tidal volumes in diseased cases versus the hypothetical healthy case. From Table 2, it can be seen that the healthy, diseased 20–40, and diseased 20–60 simulation tidal volumes were 0.592 L, 0.392 L, and 0.248 L, respectively. Thus, the results in Fig. 5 and Table 2 demonstrate that the tidal volume of the whole lung decreases as the severity of COVID-19 in the affected consolidated regions increases. It can also be seen in Table 2 that the right middle lobe shows the least difference in its tidal volume between the healthy and diseased simulations.

Fig. 6

Time-dependent volume change of the entire lung, determined at the trachea (a), and individual lobes (b–f). Values for volume change obtained by integration of flow rates.

Table 2

Tidal volume of the different lobes during one breathing cycle based on the integration of simulation results.

Lobes	Tidal Volume under different conditions (L)
	Healthy	Diseased 20-40	% change between healthy and 20-40	Diseased 20-60	% change between healthy and 20-60
Whole lung	0.592	0.392	−33.8%	0.248	−58.1%
Right upper	0.106	0.070	−33.8%	0.040	−61.9%
Right middle	0.079	0.071	−9.8%	0.060	−23.3%
Right lower	0.111	0.062	−44.2%	0.030	−72.6%
Left upper	0.196	0.137	−29.9%	0.091	−53.7%
Left lower	0.101	0.052	−48.8%	0.026	−74.0%

The volumes of GGO and consolidated regions calculated from CT images at end-inspiration were also quantified and are presented in Table 3 . As previously described, these volumes were determined directly by thresholding the CT images. Note that actual COVID-19 affected volumes are likely to be slightly larger than presented, as very small unconnected “islands” of damage had to be excluded from analysis for volume meshing purposes. In Table 3, “COVID-Afflicted %” is the combined portion of consolidated and GGO regions. It can be seen that at end-inspiration, the right middle lobe shows the lowest percentage of damage by COVID-19, followed by the left upper lobe by a large margin. The two lobes with the smallest volume of air and the highest amount of damage in both states are the left and right lower lobes (Table 3).

Table 3

Total volume of air, GGO, and consolidated regions at end-inspiration obtained from segmenting the CT image.

Volume (L)	Right Upper Lobe	Right Middle Lobe	Right Lower Lobe	Left Upper Lobe	Left Lower Lobe	TOTAL
Air	0.252	0.324	0.136	0.518	0.118	1.348
GGO	0.189	0.079	0.263	0.187	0.126	0.844
Consolidated	0.169	0.057	0.225	0.239	0.228	0.918
Total	0.61	0.46	0.624	0.944	0.472	3.110
Air %	41.3%	70.4%	21.8%	54.9%	25.0%	43.3%
GGO%	31.0%	17.2%	42.1%	19.8%	26.7%	27.1%
Consolidated %	27.7%	12.4%	36.1%	25.3%	48.3%	29.5%
COVID-Afflicted %	58.7%	29.6%	78.2%	45.1%	75.0%	56.7%

Table 4 compares the share of ventilation that goes into each lobe during tidal breathing for each of the three simulated scenarios versus the values calculated from the CT using images at end-inspiration and end-expiration. Based on these results, it can be seen that the right middle lobe which shows the lowest COVID-19 infiltration of 29.6% at end-inspiration, demonstrates a 10% increase in tidal volume percent share as the simulations progress from healthy to the diseased 20–40 and 20–60 states. The left upper lobe which has a relatively low 45.1% infiltration at end-inspiration, shows an increase in tidal volume percent share with COVID-19 progression, though not as drastically as the right middle lobe, with only a 3.4% increase (Table 4). At end-inspiration, the right upper, right lower, and left lower lobes demonstrated remarkable COVID-19-induced damage with 58.7%, 78.2%, and 75.0% affliction, respectively, and each lobe showed a decrease in tidal volume percent share with disease progression. The less-afflicted right upper lobe showed a far more minor change in tidal volume (−1.6%) than the right lower lobe (−6.5%) or the left lower lobe (−6.5%).

Table 4

Difference between the lobar share of tidal volume obtained from 4D CT (4th column), simulations of healthy and different disease states (columns 1 to 3), and Jahani et al. [48] study on lobar distribution of ventilation in healthy subjects (5th column).

	Lobar Share of Tidal Volume of Healthy simulation	Lobar Share of Tidal Volume of Diseased 20–40 simulation	Lobar Share of Tidal Volume of Diseased 20–60 simulation	Lobar Share of Tidal Volume from CT	Mean Lobar Air Volume from Jahani et al.
Right upper	17.9%	17.9%	16.3%	23.5%	20.9%
Right middle	13.3%	18.1%	24.3%	8.6%	10.7%
Right lower	18.7%	15.8%	12.2%	13.2%	22.8%
Left upper	33.1%	35.0%	36.5%	38.8%	25.8%
Left lower	17.1%	13.2%	10.6%	13.7%	20.6%

Additionally, Table 4 contains a column showing lobar air volume fraction, averaged over inspiration and expiration and converted to a percentage, as defined in a Jahani et al. [48] study investigating regional healthy lung deformation and ventilation with 4D-CT. This column was added for validation purposes to show a general agreement with our healthy model and CT-based estimations of lobar ventilation distribution. Note that lungs differ from person to person morphologically, but general trends in shape persist in most lungs; in this case, our healthy simulation successfully predicts that the left upper lobe is the largest, followed by the right lower, right upper, left lower, and right middle lobes, which mirrors Jahani et al.‘s results [48].

The pressure distributions of the lung in the three different simulation scenarios are visualized in Fig. 7 . All pressure distribution results are shown at both maximum inspiration and maximum expiration. Qualitatively, it can be discerned from Fig. 7 that the subject had the highest magnitude of pressure during both inhalation and exhalation in the right middle lobe and left upper lobe. These results correlate with the percentage of healthy (air-filled) acini determined from the CT images (Table 3). As these two lobes had the lowest percentage of COVID-19 affliction, they were expected to have the most unrestricted airflow. It is also notable that the three more COVID-19 affected lobes, the right upper and right and left lower lobes, show very different pressure distributions and magnitudes when comparing the healthy case to the diseased simulations. Furthermore, the diseased simulation scenarios show much more heterogeneity in air pressure distribution throughout the lung. Table 5 shows the average pressure at the terminal bronchioles in the entire lung for each simulation. From these values, it can be seen that the healthy lung simulation produced the highest average pressure magnitude and the least heterogeneity, while progressive disease states caused an increase in standard deviation and a substantial decrease in mean pressure.

Fig. 7

Lung pressure distribution as viewed from the front in a healthy simulation (a, b), diseased 20–40 simulation (c, d), and diseased 20–60 simulation (e, f), at maximum inhalation (a,c,e) and maximum exhalation (b,d,f). Note that the right lung appears on the left in this image and vice versa.

Table 5

Mean and standard deviation of lung pressure at the terminal bronchioles in each simulation.

Average Pressure (Pa)	Maximum Inspiration	Maximum Expiration
Healthy	−16.9 ± 1.46	17.5 ± 1.58
Diseased 20-40	−9.63 ± 2.63	8.99 ± 2.99
Diseased 20-60	−5.06 ± 3.03	4.85 ± 3.02

Discussion

Few computer modeling studies have been performed with the aim of better understanding ventilation dynamics changes in advanced cases of COVID-19. Three purely mathematical yet elegant and informative models of COVID-19 effects on pulmonary ventilation and perfusion were presented by Voutouri et al. [49], Busana et al. [50] and Herrmann et al. [51]. Another mathematical modeling study by Weaver et al. [52] simulated the effect of increased respiratory effort of patients with COVID-19 acute hypoxemic respiratory failure during spontaneous breathing on parameters associated with lung injury such as tidal swings in pleural pressure. Additionally, a computational fluid dynamics model of airflow in the upper airways of COVID-19 patients was also recently presented by Pan et al. [53]. In our simulation-based study, a virtual physics-based hypothesis-testing platform was developed and presented to study lobar ventilation dynamics of COVID-19-infected lungs. In particular, the mechanical changes in severely COVID-19-afflicted lungs compared to theoretically healthy lungs were modeled and examined in a multi-scale modeling framework. To the authors' knowledge, this is the first in silico modeling approach to use 4D CT images of COVID-19-induced lung damage to simulate regional airflow and pressure distribution in the COVID-19-infected lung. We present this study as the first step towards patient-specific physics-based models of COVID-19-afflicted lung dynamics and to lay the foundation for more detailed and individualized in silico models currently in development by our group. The multi-scale approach presented here uses patient-specific lung geometry and injury patterns obtained from CT images of a patient with advanced COVID-19 and accounts for changes in flow rate into individual lobes as a consequence of COVID-19-induced lung injury and decreased tissue compliance. Heterogeneous damage in the parenchyma of the lung due to COVID-19 was represented through changing the in surfactant levels and tissue compliance at the acinar level. In addition, our computational models enabled the visualization of acinar-level impacts of the infection on global lung dynamics. Images obtained through 4D CT were processed for efficient image segmentation and conversion to physics-based computer models. This methodology provides a solid foundation for future investigations of other potential mechanisms of COVID-19 damage to the lungs and the ensuing effects on global lung function through a computationally efficient approach.

Our disease scenario simulation results showing differences in the lobar distribution of tidal volume are in agreement with previous studies demonstrating that lung damage in advanced ARDS, which resembles Type H COVID-19 pneumonia as categorized by Gattinoni et al. [54], can decrease the compliance of injured regions to a large extent, hence decreasing airflow [55]. As expected, the least damaged lobe (right middle, 29.6% damage at end-inhalation) showed a large increase (83%) in tidal volume share when comparing the healthy and the most severe simulated disease case, and the most damaged lobe (right lower, 78.2% damage measured at end-inhalation) experienced a significant decrease in ventilation share (−35%) in the same comparison. Meanwhile, the full lung (56.7% damage at end-inhalation) also showed a notable overall reduction in tidal volume when comparing the healthy and the most severe simulated disease state (−58%).

In this study, analysis of damaged regions from the segmented CT images showed that lower lobes contained higher amounts of damage (combined GGO and consolidated regions) than other lobes (Table 2, Table 3), in this patient. These findings correspond to those from previous studies of CT imaging of COVID-19 lungs [56]. In a study using ventilation/perfusion single-photon emission computed tomography combined with computed tomography (V/Q SPECT/CT) performed in COVID-19 patients, large ventilation defects in the subpleural areas were observed [57]. The researchers observed that ventilation defects were present in the GGO areas while perfusion was largely preserved. However, in more severely damaged areas of the lung where complete alveolar filling and parenchymal lesions and fibrosis were present, perfusion defects were additionally observed which can be a sign of involvement of capillary walls [57]. Furthermore, the researchers proposed a potential adaptive mechanism where redistribution of ventilation towards the healthy parenchyma occurs [57]. This is indeed what our multi-scale model demonstrates: as the severity of damage in the right lower, left lower, and to some extent in the right upper lobe was increased, ventilation was redistributed to the least injured right middle and left upper lobes (Table 4). While our computer model in its current version only includes ventilation, we are actively developing a perfusion model which will be coupled to our ventilation dynamics model and will enable us to simulate microangiopathy induced by COVID-19.

Our simulations reasonably predicted a decrease in overall tidal volume as the level of lung damage and surfactant loss was increased. Although preserved compliance has been reported in early CARDS, lungs in more advanced COVID-19, such as those of our patient, are reported to have decreased compliance, in line with typical ARDS [54]. The decrease in tidal volume represents this reduction in compliance, as stiffer acini in the pressure range of tidal breathing would have a smaller change in volume compared to healthy acini with normal stiffness. Overall, the data in this study exhibit trends of decreased airflow to areas most affected by COVID-19 damage in an advanced case of the infection. The lower lobes showed a greater change in percent share of tidal volume in disease state simulations, while the upper lobes showed a less overall shift in percent share of tidal volumes and were less affected by COVID-19 (Table 3, Table 4). Our hypothetical simulation of healthy lung ventilation distribution among the different lobes was in qualitative agreement with the results of Jahani et al.‘s [48] 4D CT analysis of healthy lungs (Table 4). However, the ventilation distribution was notably altered in the COVID-19-infected lung simulations, with less damaged lobes receiving higher portions of tidal volume and more damaged lobes receiving a lower share of tidal volume (Table 4). Thus, regional differences seen in tidal volume distribution and subsequent changes in lobar share of tidal volume due to damage hint at how heterogenous damage in acinar regions may affect global lung function in a region-specific manner.

However, discrepancies in the lobar share of lung volumes between our CT analysis and simulations exist and are most evident when considering the right middle lobe (Table 4). The difference in values found in our study compared to 4D CT measured values could be the result of simplifications and assumptions of our model, specifically not accounting for gravity with respect to the position of the patient during imaging and the interdependence of gravity and tissue mechanics [58] or not incorporating the interdependence of acini [59] and collateral ventilation. We would like to emphasize here that the presented work is the first step toward patient-specific modeling of COVID-19 lung dynamics and hence is not intended to make quantitative predictions about ventilation dynamics, nor is it meant to serve as a clinical decision support system for clinicians at this stage. Rather our aim in this study is to lay the foundation for and take the first steps toward developing a patient-specific modeling workflow for the investigation of COVID-19-afflicted lung dynamics. We acknowledge this study has limitations that we intend to address in future studies. For example, a uniform pleural pressure was applied to all acini in our models. However, gravity has been shown to affect the spatial distribution of pleural pressure, tissue compliance, and acinar volumes [58]. Swan et al. [60] showed that even in the healthy lung, distribution of tissue compliance is spatially non-uniform. We plan on adding the effect of gravity on pleural pressure and tissue compliance distribution in future studies.

Similarly, collateral ventilation and inter-acinar interactions can affect the uniformity of pressure distributions and thus may improve the fidelity of human lung digital twins [59]. Airway deformation can also affect airflow and pressure distributions but was not included in our model; a fluid-structure interaction modeling approach would be able to capture the detailed interaction of airflow and airway wall mechanics but can be computationally more expensive compared to a reduced-order model like the one presented in our study [53,61]. Moreover, gas exchange and ventilation-perfusion coupling are important considerations for comprehensive pulmonary dynamics modeling and need to be accounted for in future studies of COVID-19-afflicted lungs [62,63]. Likewise, the lack of spirometry data of the patient and not accounting for lung motion via image registration did not allow for the implementation of patient-specific boundary conditions and direct validation of the model against 4D CT data. We are developing image registration approaches based on the work of other researchers [64,65] to account for lung motion to make the model's predictions more credible and accurate.

While surfactant dysfunction accounted for acinar level damage in this study, other forms of damage, such as diffuse alveolar damage, microangiopathy, edema, and fibrosis, lead to altered lung mechanics in the COVID-19 affected lung [36]. Including more detailed alveolar mechanics models and accounting for different types of damage in the simulations, could further elucidate the impact of the disease on global lung function. Furthermore, we modeled surfactant loss at two levels of damage based on ground-glass opacities and consolidations. Damage distribution in the lung will be more accurately represented if a continuous range of damage based on CT-derived Hounsfield values is incorporated in the model. Furthermore, we acknowledge that population-based in silico studies are bringing researchers a step closer to making virtual clinical trials a reality [66]. While the modeling workflow was demonstrated for one patient here, we aim to perform the modeling and simulation for a larger cohort of patients to study intersubject variabilities of airflow as the imaging data is collected and analyzed and our modeling workflow becomes more automated and streamlined.

In conclusion, this study presents a major step towards a modeling workflow capable of simulating the effect of heterogenous COVID-19-induced lung damage on ventilation dynamics in a patient-specific manner. The in silico model reasonably predicted redistribution of ventilation from severely damaged lung lobes to the lobes the least affected by viral insult in advanced COVID-19. This modeling study lays the foundation for patient-specific investigations of pulmonary ventilation in COVID-19 patients and individualized treatment strategies.

Conflicts of interest

The authors declare no competing interests.

Data availability

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Author contributions

A.V., E.D., S.M., and A.P., performed the method development and simulations and wrote the manuscript. M.S.P and A.J. and K.Y. collected the CT data. S.M.G. and V.M.contributed to study conceptualization and protocol development. All authors reviewed the manuscript and contributed to the scientific discussion.

Declaration of competing interest

None Declared.