A combined three-dimensional in vitro–in silico approach to modelling bubble dynamics in decompression sickness

The growth of bubbles within the body is widely believed to be the cause of decompression sickness (DCS). Dive computer algorithms that aim to prevent DCS by mathematically modelling bubble dynamics and tissue gas kinetics are challenging to validate. This is due to lack of understanding regarding the mechanism(s) leading from bubble formation to DCS. In this work, a biomimetic in vitro tissue phantom and a three-dimensional computational model, comprising a hyperelastic strain-energy density function to model tissue elasticity, were combined to investigate key areas of bubble dynamics. A sensitivity analysis indicated that the diffusion coefficient was the most influential material parameter. Comparison of computational and experimental data revealed the bubble surface's diffusion coefficient to be 30 times smaller than that in the bulk tissue and dependent on the bubble's surface area. The initial size, size distribution and proximity of bubbles within the tissue phantom were also shown to influence their subsequent dynamics highlighting the importance of modelling bubble nucleation and bubble–bubble interactions in order to develop more accurate dive algorithms.


Linear Elasticity derivation
The initial starting point for the derivation is the statement of infinitesimal strain theory as made by [1] where a derivation for σ rr (radial component of stress) can be found in Problem 2 of [1].
1 r 2 d dr r 2 ur = const = 3a (1) U σσ = a + b r 3 (4) σ rr = 3Ka − 4b mu r 3 (5) At this point boundary conditions are needed to find expressions for b and P B It is assumed that the material is linearly elastic, homogenous and that displacements are small. To solve the problem two situations must be considered, the reference configuration and a time t = t 1 later at which time the external pressure has been altered. In this derivation we are considering an infinite block of tissue which at time 0 contains a bubble with radius R 0 . As the external pressure changes the pressure at the outer boundary changes as does the stress at the bubble boundary. This causes a force imbalance at the bubble surface causing it to increase or decrease in radius until the pressure inside the bubble is once again in balance with the external pressure from the tissue. In this instance we are ignoring diffusion and just effectively looking at the boyle's law contribution to the bubble radius change as a consequence of change in external pressure.
Consider Supplementary Figure 1 Initially we have an initial external pressure of P 0 , and an initial bubble radius or R 0 . At this time there is a constant stress through the tissue of σ rr Hence at the boundary i.e. when r=∞ (2) is only valid if: and hence a = − P 0 3K So (2) becomes For when R B = R 0 At this point we use introduce the reference configuration. This is the configuration of the system if it were unstressed. To do this we say that if σ rr we zero throughout the tissue there would exist a void with a radius of R ref . The displacement U r will always be formulated as the difference between the current radius and this reference radius: U r = R B −R ref We already have an expression for this at R = R 0 give above. At later times we therefore have the displacement given as the displacement of R 0 plus the displacement between R 0 and the current R B .
From the above we also have the relation between R 0 and R ref at r = R ref as: 2.2 Situation 2: Expression for R B Consider Situation 2: If we now move to a new time point. At this time the external applied pressure P amb has changed and the stress throughout the tissue is no longer constant. We have: At r = ∞ σ rr = −P amb (14) where So we have four simultaneous equations from the above The ideal gas law gives us: to summarise; the known and unknown quantities in the above 4 equations as Now we want an expression for R B in terms of the change in external pressure P amb rearranging (22) we have Subbing (22) into (20) rearranging (21) subbing into (19) Subbing in the expression for P B and (22) Substitute the expression for b into (24) simplifying we may obtain the quartic from which R B can be found: Or by subsititution of (27) into (24) we can write the the expression in the form of the Young-Laplace equation:

Governing equations
Here we derive the governing equation for the in-silico model of bubble growth in a collagen gel tissue phantom. During the compression and time at increased pressure, the tissue phantom will saturate with dissolved gas based on Ficks law of diffusion: With the boundary conditions at the tissue edge given by Henry's Law: Where C g is the concentration of the gth dissolved gas, h g is Henry's constant and P g amb is the partial pressure of the gth gas. All gases are considered to behave as ideal gases hence: Where R is the gas constant, and T is the temperature. Alternatively where α is the specific gas constant, given by: where M is the molar mass of the gth gas. (Temp is assumed to be 37 • ). During decompression the external pressure decreases and bubbles will grow by both Boyle's Law: and by diffusion of dissolved gas into the bubble following Fick's first law: Where again Henry's Law is used as the boundary condition for (38): is the internal pressure of the bubble and this is described by the Young-Laplace equation Where Ω(R B ) is a term to account for the elastic resistance of the tissue to bubble growth. Substituting m in (38) with the perfect gas law we have d dt Expanding the LHS of (41) with the product rule and including the separate affects of both N 2 and O 2 we have d dt Assuming the mole fractions are constant throughout and using the substituting in the Young-Laplace equation to the LHS of (42) we find Where C tot is the sum of the individual gas contributions. These equations can be implemented for any number of gases. Current simulations use only Nitrogen and Oxygen. Where multiple gases are used the partial pressure of each gas is calculated from the mole fraction and the total pressure. This affects the boundary conditions and separate diffusion through the bulk is calculated for each gas.

Computational Implementation
To numerically implement the governing equations (2) and (5) the tissue phantom is defined as a 3dimensional array of points with distance between each point. Each grid point or node is described by its 3D Cartesian co-ordinate and represents a unit volume of the system. At each node the concentration of dissolved gas and phase (liquid or gas) of the node is stored. The pressure profile is discretized into time portion . At each time point the pressure external to the tissue block P amb is given by this profile. P amb is used in Henrys Law (eq. 2) to set the tissue phantom dissolved gas boundary values at each time point. From these boundary values the dissolved gas concentration at every liquid node can be calculated. This provides a 3D dissolved gas concentration, which is then used to calculate the dissolved gas gradient at the tissue phantom bubble interface. This is used in Eq. (5) to calculate the new radius of the bubble. The governing equations are discretized a solved via finite difference methods. Diffusion through the bulk of the tissue phantom is solved to second order accuracy using a forward time centered space (FTSC) scheme, utilizing a seven point stencil. The pressure change and bubble radial change are discretized to first degree accuracy,

Gas liquid interface
A system, such as the one described here, which consists of two phases, requires particular care in the treatment of the phase interface. The situation arises in many fields and there are two broad approaches: diffuse or sharp interface techniques. A diffuse interface technique assumes the phase boundary to have a finite thickness over which there is a smooth change in physical properties such as concentration. By comparison, sharp interface techniques model the boundary as a surface over which a discontinuous change in physical properties occurs. A diffuse interface technique requires a grid resolution of h i /4 where hi is the interface thickness. This makes the technique less suitable for computational domains far larger than the interface thickness. One way to overcome this restrictions is the adaptive mesh schemes to refine the grid in the vicinity of the interface, however this technique requires computationally expensive re-gridding at each time step. In our case a sharp interface technique was employed on a fixed Cartesian mesh in this model, as high computational efficiency was deemed of greatest importance given both the size of the computational domain and the time scales over which simulations were required (several hours). Due to this approach specific methods for any computations involving grid points on or crossing either the phase boundary were required. Such computations occurred in two cases: i) during bulk diffusion calculations where one of the seven stencil points fell within the bubble (see Supplementary Figure 2). ii) For calculation of the concentration gradient at the bubble boundary. In case i) the point in question is know as a boundary point, when a boundary point is reached in the calculation, the stencil point within the bubble is replaced by one at the bubble surface with concentration equal to CB. This caused the seven-point stencil to become asymmetrical and violate the stability criteria. Hence the dissolved gas concentration of the boundary point is found by inverse distance weighted interpolation from the surrounding neighbours rather than the FTSC scheme used for other bulk points. Case ii)-calculation of the gradient at the bubble surface ∂C g B ∂r | R B is dealt with using a normal probe method described in [2]. A set of spherically symmetrical points around each bubble is created as shown in Supplementary Figure 2 LHS. The concentration of these points is found using inverse distance weighting interpolation of the 6 nearest neighbours. Using the new array of spherical points a 1D forward difference approximation with 2nd order accuracy can be constructed in terms of the radial co-ordinate. where C r+2 and C r+1 are the concentrations at points r + 1 and r + 2, and r is the radial distance between points shown in Supplementary Figure 2. C B is the concentration of gas at the bubble surface as given by Henrys law, in all simulations δr = h.

Mass Conservation
Mass conservation is of particular importance for sharp interface Cartesian grid methods. Such methods are prone to local violation of mass conservation due to the movement of the boundary over grid points, creating so called fresh and dead cells which may act as spurious sources and sinks. In addition errors may arise from the discrepancy between the actual smooth curved bubble surface and its computational depiction using the cuboid volume nodes (a so-call staircase error). In order to ensure mass conservation, the mass of gas in the bubble (as calculated from the perfect gas law) is compared to the total dissolved gas concentration of the tissue and the mass flux. Total mass of gas in the tissue must be calculated accounting for the partial volume of nodes which are cut by the bubble surface as shown in Supplementary Figure 2. If a mass deficit is found, the mass is transferred to the tissue nodes directly surrounding the bubble in proportion to their volume. This method can be shown to ensure good mass conservation by simulation of an oscillatory external pressure profile with no flux boundary conditions applied to the tissue edges (see Supplementary  Figure 3).

Image Analysis
The error in bubble radius was established by measures of calibration standard metal parts, (measured with a micrometre), error in measurement of this was found to be 0.005mm. To analyse bubble growth profiles, time-lapse images of multiple bubbles within the tissue phantoms were taken starting at the commencement of decompression. The radial trajectories resulting from bubbles were fitted using a rhobust non-linear regression (Prism 6). The regression model chosen was

Parameterisation
Material parameters were set based on a literature search, separate tables of the five material parameters are shown below. In addition the parameters α N 2 and α O 2 were taken from [3]. Conversion from Ostwald coefficient (L) to Henrys constant is done via the formula from 0.014520(L) Water N 2 35 [11] 0.073(L) Sheep bone marrow N 2 37 [11] 0.015(L) Calf Brain N 2 37 [11] 0.0133(L) Olive oil O 2 37 [11] 0.0261(L) Whole man O 2 37 Table 2: Estimations of L in various biological tissues Emerson [13] pg 87. (assuming the density of biological fluid to be 1), given by where R spec is the specific gas constant for the gth gas. Values chosen for collagen gels were 0.027 for O2 and 0.0145 for N2, in the modelling literature there is also agreement with these values [14,15,16,17,18,19].

Bubble nucleation
Bubble nucleation is the subject of a separate publication however the distibution of nuclei was investigated through images of the tissue phantoms after a decompression profile by removing the phantoms from the pressure chamber and imaging them in a cross sectional view with phantoms placed in a clear petri dish containing PBS and a digital camera mounted to the side. The images were analysed manually in imageJ by taking a random rectangular slice of the whole depth of the phantom and marking the centre of all bubbles in the image section. The vertical distance of the bubble from the phantoms upper surface was normalised to the phantom thickness and plotted as a histogram and a cumulative frequency distribution shown in Supplementary Figure 5.