Article Text

## Footnotes

Contributors All authors have contributed substantially to the data generation, content and drafting of the manuscript.

Funding The authors have not declared a specific grant for this research from any funding agency in the public, commercial or not-for-profit sectors.

Competing interests None declared.

Patient consent Not required.

Provenance and peer review Not commissioned; internally peer reviewed.

## Statistics from Altmetric.com

### Key messages

Many mechanisms of injury seen in military conflict are impossible or very difficult to replicate and study in human or large animal subjects.

With comprehensive baseline data and thorough application of known physiological processes, accurate ‘in silico’ models can facilitate such research.

Such computerised models can run multiple study scenarios quickly and reproducibly and help the Ministry of Defence reduce the use of live animals in research.

Our simulator can accurately model primary blast lung injury as well as the potential effect of a range of medical interventions.

## Introduction

Recent advances in capability to accurately model and simulate human pulmonary pathophysiology have opened up the possibility of rationally ‘designing’ new multi-intervention treatment strategies in silico by exploiting the speed, reproducibility and cost-effectiveness of ‘virtual’ patient trials. In contrast to trials on both animal models and human patients, in silico models of individualised patient and disease pathology are completely configurable and reproducible—different treatments, or combinations of treatments, can be applied to the same spectrum of virtual patients, in order to understand their mode of action, quantitatively compare their effectiveness in multiple different scenarios and optimise interventions for particular clinical objectives. Such ‘virtual’ trials using simulations that are based on real patient data can offer comparable utility and validity to clinical trial data.1 2 The modelled system of interest will be, through necessity, represented in a simplified form and will inevitably require the use of reasonable assumptions in order to replace knowledge gaps.

As discussed in this paper recently,3 one such area where computational modelling could be of benefit is in improving the management of primary blast lung injury (PBLI). PBLI is part of a syndrome of injury resulting to a high explosive shock wave. It consists of a series of pathophysiological consequences resulting in a spectrum of adverse cardiovascular and respiratory consequences.4

The situation is challenging, but computational modelling could provide considerable insight into improving the understanding and treatment of PBLI, and in this paper, we provide a brief outline of our primary blast lung injury simulator, its validation and some preliminary results.

## Building a model

Model construction consists of several ordered stages. Having defined the problem or question, the important variables and properties are identified (ie, cardiorespiratory physiology in our case) and used to *build* the model. The model is then studied, validated and, finally, used. Knowledge gaps about the real-world system to be modelled will need to be filled by clearly stated assumed relationships or a ‘best guess’. The validity of the model output is limited by the validity of any assumptions made and in part helps to establish the models ‘domain of validity’. Natural laws that apply to the model should be identified and documented (box 1 for our model). Subsequently, the model must be optimised in some way such that it mirrors the system of interest as closely as mathematically possible. Common methods of mathematical optimisation include applying genetic algorithms, sequential quadratic programming and application of a mesh adaptive search. The model is then validated, whether by simple visual inspection, qualitative or quantitative means before allowing the model to become actively involved in research.

### Natural laws relevant to our model

Conservation of mass.

Universal law of gravitation.

Conservation of energy.

First-order rate laws.

Law of mass action.

Newtonian fluid laws.

## Our model

The primary blast lung injury simulator proposed here is a dynamic, discrete deterministic model built using the Matlab technical computing package (MATLAB V.8.5 R2015a). The complete simulator is an integrated model of the cardiovascular and respiratory systems, each of which has a separate signal generator completely configurable by users (Figure 1). The basic arrangement of the model has been described previously, and here we outline the modifications developed in order to create a militarily and clinically relevant simulator of primary blast lung injury. Following review of the pathophysiology of the disease process, we identified the need to add spontaneous ventilation, an enhanced and more accurate representation of lung injury constituting of alveolar rupture and pulmonary oedema formation and ventilator-induced lung injury (VILI) to the capability of our existing model. The physiological equations adopted and described are those currently accepted as valid.

## Spontaneous breathing

In order to apply our model in the prehospital environment, we needed to accommodate spontaneously breathing patients. The new spontaneous breathing module adds another signal generator which produces a sinusoidal signal to breathe which can be manipulated to replicate the known physiology of shock wave exposure (Figure 2A).

Spontaneous breathing is simulated by incorporating a variable P_{spont,} which represents the pressure generated by the respiratory muscles. The value of P_{spont} generates the force acting on each alveolar compartment (Figure 2A). P_{spont} is modelled on the profile presented by Roth *et al*.5

(1)

where is the influence of gravity and is a dynamic time-varying term. is distributed (Figure 2B) between = 3.75 and = −3.75.6 is determined through

(2)

where , , T and ∅ are parameters that can be used to generate classic spirometric breathing profiles.7 For example, a configuration given by , T=4 and can generate a sufficiently accurate representation of spontaneous tidal breathing simulation (Figure 2B).

## Ventilation perfusion mismatch caused by blast lung injury

The fundamental physiological disturbance leading to impaired gas exchange following shock wave exposure is ventilation-perfusion (V/Q) mismatch. Therefore, a key requirement for the simulation of pulmonary pathology including PBLI is the ability to accurately reflect and model the V/Q mismatch observed in our source data. In our solution, the V/Q mismatch can be manipulated by simulating recruitment and derecruitment of alveolar compartments and controlling the rate of perfusion across each of the individual alveolar compartments in the model.

The total pulmonary vascular resistance (PVR) is determined by

(3)

where is the number of alveolar compartments (set to 100), and the resistance for each compartment is defined as

(4)

is the default vascular resistance for the compartment with a value of 160 N_{A} dynes.s.cm^{−5}.min^{−1} (all resistances are in parallel, giving a default total pulmonary resistance of 160 dyn.s.cm^{−5}.min^{−17}).
is the vascular resistance coefficient used to implement the effect of hypoxic pulmonary vasoconstriction.

The model includes the effect of radial compressive and axial stretching forces exerted onto pulmonary capillaries as a result of increase in lung volume and pressure. The overall effect on resistance to flow through each capillary is difficult to quantify, but we assume the following: (1) at alveolar volumes above the functional residual capacity (FRC), the vessels become compressed and raise the PVR, (2) at alveolar volumes below FRC, the vessels can collapse and thus result in an increase in PVR, while closer to FRC the PVR remains unaffected. The resultant ‘U’ shape change in PVR at around the FRC has been suggested previously and has been implemented in this model as follows. The PVR is determined as given in equation (4), but the vascular resistance for each alveolar compartment, , has been modified to

(5)

is calculated as follows:

(6)

where, is the pressure generated within the i th alveolar compartment, is the volume of the i th alveolar compartment, is a constant representing the volume of the alveolar compartment at rest (fixed to 3000/ equating to a functional residual capacity of 3000 ml). and is used to adjust the effect on pulmonary vascular resistance. Based on earlier studies,8 is set to 1 and has been set to 30 but they can be modified to fit patient data.

In addition to the effect on PVR, the average alveolar compartment pressure within the lung exerts an extrinsic pressure which is applied to the intrathoracic vascular compartments. This phenomenon is known as splinting. The pressure calculation of the compartments within the thoracic cavity therefore has an additional term, , added to them, representing the intrathoracic pressure

(7)

for extrathoracic compartments. Within the thoracic cavity, a range of values (0.1–0.8) is used for to fit patient data if available.

One of the potential complications of mechanical ventilation is VILI. VILI is characterised by serious lung parenchymal injury including increased permeability of the alveolar-capillary membrane with resultant pulmonary oedema and surfactant inactivation. This leads to atelectasis and systemic biotrauma. VILI is associated with overstretching of alveolar units to volumes well above the resting volume, cyclical opening and closing of alveoli and pulmonary epithelial injury due to barotrauma. Acutely injured lungs are particularly prone to developing VALI. The simulator can be used to study strain and compliance at the alveolar level caused by mechanical ventilation with all of the clinically relevant the tidal volumes, pressures and modes of ventilation. The most important components of this part of the model are the alveolar threshold opening pressure (TOP)9 and the time required for this to be achieved (alveolar recruitment). Recruitment in injured lung is a time-dependent process, with different airways recruiting at different times, once the threshold pressure has been achieved. As usual in this model, the equations modelling this are solved iteratively as a discretised system in 10 ms intervals. This time-dependent recruitment phenomenon is achieved by the introduction of a parameter
. For collapsed compartments,
is set to
which represents the time it could take for collapsed alveoli to open after a threshold pressure is reached. Once
is satisfied, the counter
decrements during every iteration and triggers the opening of the airway (
= 1) as
. Otherwise
is set to a high value (10^{10}) to represent a collapsed airway. We based the range of values for TOP used in these simulations on the work done by Crotti and collaborators10

## Pulmonary oedema formation

An increase in extravascular lung water is a common occurrence to the pathophysiology of blast lung injury and other forms of acute lung injury. Our most recent addition to the simulator is the modelling and validation of a mechanism to simulate the production and effect of focal or pericontusion pulmonary oedema. In health, interstitial fluid volume is maintained by a balanced fluid exchange between capillaries and the lymphatic system. Fluid diffuses into the interstitial space by filtering through the capillaries. This fluid is removed from the interstitium into the lymphatic system via a well-known filling and pumping mechanism. In health, the rate at which fluid is added to the interstitium from the capillaries is equal to the rate at which fluid is removed into the lymphatic system. The flow rate is governed by the Starling-Landis hypothesis, under which, the fluid movement due to filtration across the wall of a capillary is dependent on the balance between the hydrostatic pressure gradient and the oncotic pressure gradient across the capillary. The rate of fluid removed by the lymphatic system is calculated in the model by the following mathematical formulae describing lymphatic flow.

If denotes the rate at which fluid is added to the interstitium from the capillaries and is the rate at which fluid is removed from the interstitium into the lymphatic system. In a normal healthy individual, . Oedema arises when the influx of fluid into the interstitium exceeds outflow ( ), producing an accumulation of interstitial fluid.11 This extra fluid results in an increase in interstitial fluid pressure and volume.

We model the change in volume ( ) of the interstitial fluid with the time slicing approach,

(8)

where N is the number of alveolar compartments, is the initial interstitial volume before the iteration starts and h is the size of the time period in this iteration. The rate is calculated by the following formula based on the Starling-Landis hypothesis which states that the fluid movement due to filtration across the wall of a capillary is dependent on the balance between the hydrostatic pressure gradient and the oncotic pressure gradient across the capillary. Using the Starling-Landis hypothesis, we calculate the rate, , via the following mathematical equation:

(9)

where is the microvascular filtration coefficient, is the hydrostatic pressure of fluid in the pulmonary capillaries and is the interstitial pressure (ie, the hydrostatic pressure of fluid outside the capillary). and are, respectively, the capillary oncotic pressure (due to the presence of non-permeating solutes inside the capillaries) and interstitial oncotic pressures (due to the presence of non-permeating solutes outside the capillaries). indicates the permeability of the blood capillaries to non-permeating solutes (such as plasma proteins).

(10)

The alveolar pressure and volume relation is represented by equation (11).

(11)

The constant
controls the intra-alveolar pressure for a given volume.
indicates the ‘constant collapsing volume’ at which an alveolus is totally collapsed. The extra fluid due to oedema results excess pressure in the alveoli which can be modelled by modifying the equation (11) as below adding a fluid compartment V_{ext},

(12)

The new term in equation (12) (V_{ext,i}) reflects increasing fluid in the alveoli and the subsequent increase in TOP and thus atelectasis. This in turn causes an increase in PVR.

## Determining the parameters of the model to create virtual patients

In order to determine the values of model parameters so that its responses match our source data, we applied the following two processes. First, a cost function that establishes the accuracy of the model response for a given parameter configuration and second, a mechanism to intelligently determine the best parameter configuration, that is, a specific parameter that gives the most accurate model responses in comparison to the source data. For the first process, we used a standard weighted and aggregated cost function.12–14 For the second process, we chose to apply a genetic algorithm (GA) provided within the Matlab Global optimisation toolbox. GAs are ideally suited, due to their ease of application and consistent performance in problems where the number of parameters could be large and consequently the parameter space could be complex. With their probabilistic and parallel nature, GAs are generally capable of converging to the global optimum solution (in this case, the best model configuration that would generate the closest possible match of the model responses to data) even in highly complex parameter spaces.

GAs are based on evolutionary concepts such as selection, mutation, recombination, etc. A randomly selected population of candidates (first generation) undergoes a repetitive process of reproduction, where selection is based on the value of the objective function (also called the fitness). Every generation, the best candidates from the previous generation (elitism) and candidates obtained through mutation and crossover, recombine to form a new population. The average fitness of the individuals in the population is expected to increase as strong individuals are protected and combined with one other and weak individuals are discarded. To speed up the parameter identification process, a parallelised computer code implementation of a genetic algorithm was employed with accelerated hugely by distributing the tasks to multiprocessors (multiple cores and/or multiple machines). Initial model calibration and analysis were performed on a 64-bit Intel Core i7 3.7 GHz PC, running Matlab (R2014a). Model calibration to data was performed using the ‘Minerva’ high performance computing cluster provided by the University of Warwick (396 nodes, each with 2× hexa core 2.66 GHz 24 GB RAM) running Matlab (2015a) with global optimisation and parallel computing toolboxes.

## Comparison data and preliminary model results

Once optimisation was complete, we exposed six virtual patients to the same blast insult as our source data and then compared averaged model predicted outcomes with actual outcomes recorded over an 8-hour period. An example of this comparison is demonstrated in Figure 3 in which values for partial pressure of arterial oxygen (PaO_{2}) and partial pressure of arterial carbondioxide (PaCO_{2}) can be seen. Our model’s predicted outcome is plotted in blue. All of the variables measured demonstrated a similar degree of agreement. The actual values are detailed in Table 1.

We are currently in the early stages of applying our model to a variety of clinical situations of relevance to the military population. Optimal mode of mechanical ventilation is of particular interest to the deployed anaesthetic and intensive care cadre and we are beginning to achieve results in a study comparing conventional low-tidal ventilation with a basic airway pressure release mode of ventilation. Preliminary results for PaO_{2} and PaCO_{2} following 8 hours of ventilation starting 1 hour after injury can be seen in Figure 4. As can be seen, our model suggests that airway pressure release ventilation results in early improvements in gas exchange with improved oxygenation and somewhat surprisingly improved ventilation. Following further refinement, we will extend this ventilation study to a 48-hour interval period.

## Limitations of the model

Our model assumes prior good health and sea level atmospheric conditions. Autonomic reflexes are neglected because, in the studies used for model calibration, it is assumed that the cardiovascular side effects of the drugs used for sedation suppressed normal cardiovascular system baroreceptor reflexes (these studies consistently reported no significant changes in heart rate throughout their interventions).15–17 We do not start modelling activity until 1 hour after injury at which time the only remaining pertinent autonomic blast response is reduced systemic vascular resistance which we account for. Effects due to increased cytokine presence in the systemic circulation (biotrauma) due to alveolar-capillary membrane damage are not currently included.

## Conclusion

Modelling an acute injury and consequent disease process is a complicated process that requires a thorough knowledge of the relevant physiology and pathophysiology. We have described here the core components of our blast lung injury simulator, how it was generated and delineated some of the assumptions and natural laws we have applied to create it. The model is now in use and hopefully will be a valuable clinical tool to clinicians both at home and abroad.

## Footnotes

Contributors All authors have contributed substantially to the data generation, content and drafting of the manuscript.

Funding The authors have not declared a specific grant for this research from any funding agency in the public, commercial or not-for-profit sectors.

Competing interests None declared.

Patient consent Not required.

Provenance and peer review Not commissioned; internally peer reviewed.

## Request Permissions

If you wish to reuse any or all of this article please use the link below which will take you to the Copyright Clearance Center’s RightsLink service. You will be able to get a quick price and instant permission to reuse the content in many different ways.