Abstract
The transport of long gas bubbles, suspended in liquid, through symmetric bifurcations, is investigated experimentally and theoretically as a model of cardiovascular gas bubble transport in air embolism and gas embolotherapy. The relevant dimensionless parameters in the models match the corresponding values for arteries and arterioles. The effects of roll angle (the angle the plane of the bifurcation makes with the horizontal), capillary number (a dimensionless indicator of flow), and bubble volume (or length) on the splitting of bubbles as they pass through the bifurcation are examined. Splitting is observed to be more homogenous at higher capillary numbers and lower roll angles. It is shown that, at nonzero roll angles, there is a critical value of the capillary number below which the bubbles do not split and are transported entirely into the upper branch. The value of the critical capillary number increases with roll angle and parent tube diameter. A unique bubble motion is observed at the critical capillary number and for slightly slower flows: the bubble begins to split, the meniscus in the lower branch then moves backward, and finally the entire bubble enters the upper branch. These findings suggest that, in large vessels, emboli tend to be transported upward unless flow is unusually strong but that a more homogeneous distribution of emboli occurs in smaller vessels. This corresponds to previous observations that air emboli tend to lodge in the upper regions of the lungs and suggests that relatively uniform infarction of tumors by gas embolotherapy may be possible.
 gas embolotherapy
 air embolism
 tumor infarction
 perfluorocarbon
 microbubble
the motion of a gas bubble through a bifurcation is of particular interest in understanding air embolism and in designing gas embolotherapy strategies. Air embolism can have dangerous, even fatal, effects and can be caused by clinical procedures if the appropriate measures are not taken. The formation of bubbles in the systemic circulation is due to the bronchopulmonary pressure gradient that allows the transmission of air to the pulmonary veins. The bubble then travels through the vasculature, and, if not treated correctly, it may become lodged in the microcirculation (1, 40), causing ischemia. This condition can be fatal, if an air bubble lodges in the arterioles or capillaries in the brain or in the coronary circulation.
A number of researchers have investigated the distribution of air emboli in the lungs using animals and benchtop models (12, 13, 41). They concluded that the upper regions of the lung receive the majority of bubbles due to the buoyancy of bubbles in the blood (12). Souders (41) demonstrated that blood flow may affect the distribution of bubbles, although they did not specify the relationship between flow and distribution. A model study performed with bubbles whose diameters were much smaller than the tube diameter demonstrated that the position of the bubbles in the parent tube of the bifurcation determined which daughter branch the bubble entered (13). The majority of the bubbles were in the top onehalf of the parent tube because of buoyancy, causing the upper branch of the bifurcation to receive more bubbles than the lower branch did. The present investigation considers the transport of a single, long bubble passing through a bifurcation and provides an explanation of some of the bubble distribution behavior noted in previous studies (12, 13, 41).
Our interest in cardiovascular bubble transport is also motivated by gas embolotherapy. Embolotherapy involves the occlusion of blood flow by foreign objects for therapeutic reasons, primarily to starve tumors (2, 18, 36, 38). Previous embolotherapy methods have primarily used solid emboli and involve complicated procedures, requiring either surgical exposure of vessels near the target region or very selective catheter placement to minimize embolization to normal tissue. Currently, embolotherapy is used as a last resort after conventional clinical methods to eradicate tumors (chemotherapy, radiotherapy, surgery, etc.) have failed. Hepatocellular carcinoma and renal carcinoma are well suited for treatment by embolotherapy because surgical resection is difficult and chemotherapy and radiation treatments are often not successful (38).
Our laboratory is developing a novel type of gas embolotherapy (6, 34, 35, 44) that is less invasive and allows more precise delivery of emboli than conventional methods. This embolotherapy approach utilizes perfluorocarbon (PFC) gas bubbles that originate as liquid droplets that are small enough to pass through capillaries. When the droplets reach the target region, they are vaporized using ultrasound. Droplet preparation and acoustic activation are described elsewhere (6, 34, 35, 44). The resulting bubbles are long compared with the vessel diameter and, with volumes 125–150 times the droplet volume, are large enough to lodge in the circulation. This technique has the distinct advantages that the dose of PFC droplets can be injected in a convenient location rather than only near the tumor, and that the droplets can then be vaporized in close proximity to the tumor. In other studies, our laboratory has investigated the potential for this vaporization process to injure or rupture vessels using computational modeling (44). It is expected that the extended occlusion provided by the PFC bubbles or repeated application may be sufficient for tumor necrosis.
In this study, we investigate the splitting of long bubbles in model artery and arteriole bifurcations. Previous work has considered small bubbles and droplets traveling through bifurcations (13, 37), but none has examined the effects of gravity on the splitting of long bubbles in a bifurcation. The splitting of liquid plugs, rather than gas bubbles, in bifurcations without gravitational variations has been examined in relation to surfactant delivery to the lungs (8), and the transport of these liquid plugs has been studied in vivo on the organ size scale for pulmonary surfactant and liquid PFC delivery (7, 9). The effects of gravity on bubble transport have been investigated in straight tubes in the context of air embolism (10, 11, 17). The work presented here elucidates the relative effects of flow and gravity on bubble splitting in a single bifurcation. Bubbles may travel through several generations of bifurcations before sticking, and the splitting behavior will determine the homogeneity of flow occlusion. Previous studies have suggested that bubble sticking occurs in vessels smaller than those modeled here (40), and, therefore, we do not consider sticking in this work. The dynamics of sticking and the pressure difference at which bubbles stick are the topics of other investigations in our laboratory. Others have investigated adhesion of air emboli in the microcirculation and demonstrated the potential for surfactants and PFCs to dislodge emboli (3, 10, 19–21, 42, 43). Understanding the splitting behavior of bubbles, which is the focus of this paper, will be essential for developing gas embolotherapy strategies and provides insight regarding the optimal location for acoustic droplet vaporization (34) to achieve homogenous embolization of tumors.
MATERIALS AND METHODS
The experimental setup is shown in Fig. 1. Two bifurcations were used: an artery and an arteriole model. The daughter tubes of each bifurcation were symmetric and formed a Y shape. The size of each model was chosen to maintain appropriate nondimensional parameters for either arteries or arterioles (see Table 1). A syringe pump (Harvard Apparatus, PHD 2000) was connected to the parent tube of the bifurcation and provided a constant flow rate (Q) rather than an imposed driving pressure. This facilitated investigation of the competition between Q and gravity at the bifurcation level and ensured that bubbles did not stick. The two daughter branches were maintained at the same pressure by connecting them to two different reservoirs with the same hydrostatic pressure. The hydrostatic pressures were maintained at a constant value by completely filling the reservoirs and allowing the water to run over into a large container. The water level in the reservoirs was higher (6.35 cm) than the bifurcation. The reservoirs were connected to the bifurcation with two different kinds of tubing (Manifold Pump tubing inner diameter, 0.01 in., and clear polyvinyl chloride tubing, 3/16 in.) and a connector between them. To ensure that no bubbles were trapped in the connectors, the connectors were flushed with water (by running the syringe pump at a very high Q) before each experiment.
While the core fluid (water) was stationary, a long bubble of air, several times the diameter of the parent tube, was injected with a syringe into the parent tube before the bifurcation. After the bubble was placed in the tube, the syringe pump was turned on, and a constant Q of water pushed the bubble through the bifurcation. The tubing length before the bifurcation was sufficiently long to establish that the bubble was in steady flow before it reached the bifurcation. A digital camera (Sony, DCRPC110) acquired video images of the bubble as it moved through the bifurcation. The bifurcation was placed on a platform that could be adjusted to different roll angles. The camera was always placed at a 90° angle with the bifurcation using a Tsquare to ensure accurate alignment. Additionally, the image was checked for distortion that would occur if the camera were not properly aligned with the bifurcation, e.g., one region of the image would appear in focus, but other regions would not, if the camera were misaligned.
The video was transferred to a Pentium 4 personal comupter (Dell) via Adobe Premier and an IEEE 1394 firewire connector and converted to an ordered sequence of bitmap images for analysis. From these images, the lengths of the bubbles in the daughter branches were measured, and then the ratio of lower branch length (L_{b}) to high branch length (L_{t}) (splitting ratio = L_{b}/L_{t}) was calculated. This process was repeated in a randomized manner for the range of parameter values to obtain five measurements of splitting ratio for each setting. The data for the same parameter values were averaged, and the standard error of the mean was calculated.
The following relevant, dimensionless parameters were calculated for the experiments to allow comparison of the experiments with the physiological setting.

1) The capillary number, Ca = Uμ/σ, indicates the ratio of viscous force to surface tension force, where U is the average velocity of the core fluid (water) in the parent tube, μ is the viscosity of the core fluid, and σ is the surface tension between air and water.

2) The Reynolds number, Re = UD_{p}ρ/μ, indicates the ratio of inertial force to viscous force, where D_{p} is the diameter of the parent tube of the bifurcation, and ρ is the density of the core fluid (water).

3) The Bond number, Bo = ΔρgD/σ, represents the ratio of gravity force to surface tension force, where Δρ is the difference between the core fluid density and the bubble density, g is the acceleration due to gravity, and σ is the surface tension between air and water. Defined in this classical manner, Bo does not account for the effect of varying roll angle θ on splitting. Because the component of gravity that acts in the plane of the bifurcation is g·sin(θ), we also consider a modified Bo, Bo* = Bo·sin(θ), to account for θ. When θ = 0, we do not expect gravity to affect bubble splitting, i.e., Bo* = 0. Although Bo in these experiments is higher than in actual arterioles, Bo* is the parameter that must be matched to physiological values in order for this model system to apply to arterioles. By using small values of θ, we can mimic the range of Bo* found in human physiology.

4) Dimensionless bubble length, L_{d}, is the dimensional bubble length divided by the parent tube diameter.
The Q was varied to vary the values of Ca and Re. Changing Q changes both Ca and Re, as the mean velocity of the core fluid in the parent tube appears in the numerator of both of these parameters.
Theoretical model.
To better understand the mechanisms responsible for the bubble behavior observed in the experiments, we developed a theoretical model whose results we compared with the experimental results. The model geometry is shown in Fig. 2. As depicted, the bubble is splitting with part of the bubble entering the upper daughter tube, denoted by 1, and part of the bubble entering the lower daughter tube, denoted by 2. P_{1} and P_{2} are the pressures in the liquid (water) at the nose of the bubble in the upper and lower tubes, respectively. The velocities of the bubble noses in the upper and lower daughter tubes are denoted by V_{1} and V_{2}, respectively. The elevations of the bubble noses in the two tubes are denoted by z_{1} and z_{2}, respectively. These can be related to the distance, l_{1} or l_{2}, of the bubble nose from the end of the daughter tube and θ by geometrical relations. The bubble is considered to be close enough to the parent tube wall that liquid flow between the bubble and wall is negligible at the low capillary numbers considered here (5, 27, 28, 30, 31). Assuming quasisteady motion of the bubble, the pressure difference between the bubble and the surface of either reservoir is due to 1) the pressure jump across the bubble interface, ΔP, 2) elevation differences, and 3) viscous losses. We can express the bubble pressure, P_{b}, in terms of these contributions along either the upper or lower daughter branch to obtain the following two equations. (1) (2) where subscripts 1 and 2 denote points 1 or 2, e.g., upper or lower daughter tube, respectively; reservoir surface is denoted by the subscript e for exit; γ = ρg is liquidspecific weight; P is pressure; and h_{l} is head loss due to viscous effects in the tubes. This pressure difference calculation is similar to the analysis used by Ghadiali and Gaver (27) for the steady motion of semiinfinite bubble in a single liquidfilled tube. Here, we have neglected inertial effects because they are small compared with the effects of viscosity, surface tension, and gravity. The Poiseuille flow in the majority of liquid ahead of the bubble is unidirectional and, therefore, does not contribute to convective inertia. Note that the P_{b} is assumed uniform throughout the bubble, but varies with time as the Q is maintained.
The h_{l} between the nose of the bubble and the reservoir, for either the upper or lower daughter tube, was calculated based on viscous losses in the daughter tube and in the subsequent Tygon tubing connecting the bifurcation with the reservoir. Considering laminar flow and neglecting minor losses due to bends, fittings, and the entrance to the reservoir, since they are small compared with the viscous losses over the length of the tubes, the h_{l} values in the upper and lower daughter tubes were estimated by (3) (4) respectively, where d is the diameter of the daughter tube, D is the diameter of the Tygon tubing, L is the length of the Tygon tubing, and V is the speed of the bubble nose. The distance l varies with time, decreasing as the bubble propagates through the bifurcation. The subscripts 1 and 2 denote the upper and lower daughter tubes, respectively. The V of the bubble nose indicates the rate at which the length l decreases, so that (5) (6) The pressures at points 1 and 2, P_{1} = P_{b} − ΔP_{1} and P_{2} = P_{b} − ΔP_{2}, respectively, are related to P_{b} by the normal component of the interfacial stress boundary condition. We approximate this by the expression determined by Bretherton (5) for the pressure jump across the interface of a semiinfinite bubble in the limit of small Ca (7) where i denotes either 1 or 2, and V is the velocity of the bubble tip in the daughter tube.
Subtracting Eq. 2 from Eq. 1 yields (8) Note that the elevation difference between the bubble noses in the upper and lower daughter branches can be determined by geometry to be (9) where l_{1}(t = 0) = l_{2} (t = 0), the 39° appears due to the bifurcation angle of 78°, and it is assumed that the bubble nose does not reach the horizontal section of the daughter tube (see Fig. 2) before splitting occurs. Mass conservation requires that (10) To see the relative importance of gravity and flow, represented by Bo and Ca, we nondimensionalized the governing equations by the following scales. (11) where t is time, and dimensionless variables are indicated by circumflexes. Substituting Eq. 7 for P_{1} and P_{2} into Eq. 8 and nondimensionalizing yields (12) (13) Note from Eq. 9 that the Bo term in Eq. 12 depends on the effective Bo, Bo* = Bo·sinθ. Equations 12 and 13, along with the dimensionless versions of Eqs. 5 and 6, form a system of equations for V̂_{1}, V̂_{2}, l̂_{1}, and l̂_{2}. These can be solved by differentiating Eq. 12 with respect to time to obtain an ordinary differential equation for velocity and then solving the system of differential equations using the finite difference method in time. We solved these equations for each set of parameter values used in the artery and arteriole models. The values of l̂_{1} and l̂_{2} were determined at each time step.
To determine the splitting ratio, we determined the values of l̂_{1} and l̂_{2} at the time, t_{f}, corresponding to when a bubble of prescribed volume splits based on the imposed Q. (14) Equation 14 becomes t̂_{f} = L_{d} when nondimensionalized. The theoretical splitting ratio was then calculated based on the resulting bubble lengths at time t̂_{f}. (15) We compared this theoretical splitting ratio with the experimental splitting ratio at the same Ca and θ. Additionally, we investigated the timedependent splitting behavior by plotting L̂_{b}(t̂) and L̂_{t}(t̂) vs. t̂. Note that, in this notation, L_{b}/L_{t} without t in parentheses and without circumflexes denotes the splitting ratio at the final time, at which the bubble has split. We explicitly denote the timedependent L_{d} in the upper and lower daughter tubes as L̂_{b}(t̂) and L̂_{t}(t̂), where L̂_{b}(t̂) = l̂_{2} (0) − l̂_{2} (t̂) and L̂_{t}(t̂) = l̂_{1} (0) − l̂_{1} (t̂).
Parameters for artery model.
The diameters of the artery model bifurcation were 0.5 cm for the parent tube and 0.34 cm for the daughter tubes. These values were based on the data of Huang et al. (33), who measured the morphometry of human pulmonary vasculature. The angle of the bifurcation is 78°, approximately the angle of human pulmonary arterial bifurcations (13, 32). The ranges of the nondimensional numbers for the artery bifurcation are as follows: Ca = 1.09 to 8.76 × 10^{−4}, Re = 44.2–354, and Bo = 8.95 × 10^{−1}. Bo* for this model follows the equation Bo* = Bo·sin(θ), and the values for Bo* range from 0 at θ = 0° to 8.95 × 10^{−1} at θ = 90°. Re and Ca correspond to Q, ranging from 10 to 80 ml/min, in increments of 10 ml/min, and can be calculated by multiplying Q by a constant: Re = 4.42Q and Ca = 1.09 × 10^{−5}Q. The maximum Ca in these experiments was limited by the time resolution of the video camera. Splitting ratios were measured at θ ranging from 0 to 90°, for air bubbles of L_{d} ∼2. Typical parameter values based on the human vasculature (14, 29) and the properties of PFC gas (34, 35) are listed in Table 1.
Parameters for arteriole model.
The diameter of the parent tube for the arteriole model was 0.1 cm, with a daughter tube diametertoparent tube diameter ratio of ∼0.78, which is approximately that of human arterioles (25). The branching angle is the same as that of the artery model (78°). An experimental protocol similar to the one described for the artery model was implemented with L_{d} = 20, because bubbles are expected to be much longer in arterioles than in arteries. To match dimensionless parameters to those in human arterioles, the Q values were smaller than Q values for the artery model and ranged from 0.04 to 0.3 ml/min. This corresponds to the following ranges of dimensionless numbers: Ca = 1.09 to 8.21 × 10^{−5}, Re = 8.85 to 6.63 × 10^{−1}, and Bo = 3.58 × 10^{−2}. Ca and Re can be calculated in terms of Q by Re = 22.11 Q and Ca = 2.736 × 10^{−4} Q. The values of the dimensionless parameters in the model are similar to those of arterioles and small arteries in the human body (see Table 1). The value of Bo in the arteriole model was higher than Bo in the actual arterioles. The importance of gravity relative to surface tension depends on the θ and is indicated by Bo* = Bo·sin(θ). Thus we are able to consider physiologically relevant values of Bo*, even though Bo is different from the physiological value. Bo* ranged from 0 at θ = 0° to 1.8 × 10^{−2} at θ = 30°. The maximum value of θ used in the arteriole model was 30° because the bubble did not split at higher values, e.g., the splitting ratio was zero for θ > 30°.
To assess the effect of bubble length on the splitting ratio, we performed similar experiments with three different values of L_{d}: 6, 12, and 20. Five experiments were conducted for each L_{d}, and the results were averaged. To ensure the accurate values of initial bubble length, marks were placed on the surface of the bifurcation that corresponded to the three lengths. An individual bubble was injected until it reached the desired mark. The effect of bubble length was only considered in the arteriole model because it was difficult to place a larger bubble in the parent tube of the artery model. Additionally, it is not expected that air emboli or the bubbles created in gas embolotherapy will be larger than the bubbles used in the artery model.
RESULTS
Artery model.
Because the liquid Q was constant, the bubbles in the artery model did not stick in the bifurcation. The bubbles either split or went entirely into the upper branch. At θ = 0°, the Bo* was 0, and the bubble (on average) split evenly. These results confirmed that the pressure difference between the two daughter tube exits was zero, as even splitting would not be expected otherwise. The experimental and theoretical values of the splitting ratio L_{b}/L_{t} vs. Ca are plotted for a variety of θ values in Fig. 3. The symbol indicates the average of the five experiments at each condition, and the error bars indicate the SE. As shown in Fig. 3, there is a critical capillary number, Ca_{crit}, for all θ > 0, in both the experiments as well as in the theoretical model. When Ca < Ca_{crit}, the bubble does not split and travels entirely into the upper tube. The splitting ratio was zero in each of the individual experiments for Ca < Ca_{crit}, and, consequently, the SE was zero. Increasing θ decreases L_{b}/L_{t} at values of Ca > Ca_{crit}. For Ca > Ca_{crit}, the bubble split, and increasing θ decreases L_{b}/L_{t} at a given value of Ca. Ca_{crit} increases with increasing θ and reaches a value of 4.38 × 10^{−4} for the experiments and a value of 4.89 × 10^{−4} at θ = 90° (Bo* = 8.95 × 10^{−1}). At a given value of θ > 0, increasing Ca for Ca > Ca_{crit} increases L_{b}/L_{t}. As Ca increased further, the slope of L_{b}/L_{t} vs. Ca was not as high. The agreement between the experiments and theory was relatively close.
An interesting and unexpected bubble behavior in the experiments was noted for θ > 30°. As the bubble entered the bifurcation, it began to split. After the bubble entered the bifurcation, but before the rear meniscus of the bubble had reached the bifurcation, the lower meniscus moved backward, and the bubble went entirely into the upper tube. An example of this phenomenon is shown in Fig. 4 for θ = 30° and Ca = 1.09 × 10^{−4}. In Fig. 4A (t = 0), the bubble was transported toward the bifurcation. The bubble then began to split when it reached the bifurcation, as shown in Fig. 4B (t = 1 s). Later, more bubble volume went into the upper branch, but both menisci continued forward (see Fig. 4C, t = 1.8 s). The meniscus in the lower branch then moved backward toward the bifurcation (see Fig. 4D, t = 2.5 s). The entire bubble finally went into the upper branch as shown in Fig. 4E (t = 3.4 s). This behavior did not occur for Ca > 3.28 × 10^{−4}, but it happened more frequently as θ increased. This behavior was observed at values close to Ca_{crit} for each θ ≠ 0. On some occasions, the lower meniscus reversed direction and began to return to the bifurcation, but, as it moved backward, the bubble snapped off, the majority of the bubble went into the upper branch, and a small bubble traveled through the lower branch. This occurred infrequently in the artery model and was only at the Ca immediately greater than Ca_{crit}. The smaller bubble that snapped off did not contact the tube wall.
To examine the bubble reversal phenomenon with the theoretical model, we plotted L̂_{b}(t̂) and L̂_{t}(t̂) vs. t̂ for four Ca values close to the theoretical value of Ca_{crit} at θ = 45°, as shown in Fig. 5. For Ca = 3.33 × 10^{−4}, which is below Ca_{crit}, the bubble initially entered both branches, and L̂_{b}(t̂) increased until it reaching a maximum near t̂ = 0.7. Then L̂_{b}(t̂) decreased with time until reaching zero, indicating that the entire bubble went into the upper daughter branch and resulting in a splitting ratio of zero. For Ca = 4.68 × 10^{−4} > Ca_{crit}, L̂_{b}(t̂) also reached a maximum and began decreasing, indicating the return of the bubble meniscus in the lower branch to the parent tube. However, the bubble split before L̂_{b}(t̂) reached zero, resulting in a nonzero splitting ratio when Ca > Ca_{crit}. When Ca = 9.44 × 10^{−4}, which is much higher than Ca_{crit}, L̂_{b}(t̂) did not reach a maximum value before the bubble splits, indicating there was no bubble reversal at this Ca.
Arteriole model.
The results were similar to those of the artery model. L_{b}/L_{t} increased with Ca and decreased with θ, as shown in Fig. 6. Compared with the artery model, L_{b}/L_{t} for a fixed value of θ did not change as much with Ca over the range of Ca investigated. At θ = 30°, there was a Ca_{crit} below which the bubble did not split. However, for the other values of θ (<30°), we did not observe a Ca_{crit} and the bubble split for the entire range of Ca we investigated. Values of θ > 30° were not investigated because the bubble did not split and went entirely into the upper tube. The theoretical model overestimated L_{b}/L_{t} for each value of θ. The theoretical model results indicate a Ca_{crit} for each value of θ.
Unlike in the bubbles in the model artery experiments, the bubbles in the arteriole experimental model experiments sometimes exhibited stickslip behavior, in which one meniscus alternated between moving and sticking. As the bubble entered the bifurcation, it started to divide. Then the meniscus in the lower branch appeared to stop, and the meniscus in the upper branch increased in speed when this happened, as would be expected from conservation of mass. During this process, the rear meniscus in the parent tube continued to advance steadily due to constant Q. When the rear meniscus reached the bifurcation, the bubble split into two bubbles, and the bubble in the lower branch resumed forward motion. The bubble reversal phenomenon described above for the artery model in which the bubble entered the lower branch and subsequently returned to move entirely or partially into the upper branch was also observed in the arteriole model, but only at θ = 30° and only at Ca_{crit}. Similar behavior with the additional feature that the bubble snapped off, e.g., fragmented into two bubbles, near the bifurcation while the meniscus in the lower branch was moving backward was observed at the first Ca higher than Ca_{crit}. The portion of the bubble in the lower tube when snapoff occurred was subsequently transported through the lower branch.
Figure 7 shows L_{b}/L_{t} vs. Ca for three different average values of L_{d} (6.26, 12.71, and 21.00) at θ = 15°. For this range of L_{d}, the experiments showed that L_{b}/L_{t} decreased with L_{d}, for Ca ≥ 2.19 × 10^{−5}. As L_{d} was increased, more of the bubble went into the upper branch, but in the lower branch there was less of an increase. At each Ca value, the splitting ratio for L_{d} = 6.26 was higher than that of L_{d} = 12.71 and L_{d} = 21.00. The theoretical results exhibit similar behavior. Near Ca = 2.19 × 10^{−5}, there is a crossover of the L_{b}/L_{t} vs. Ca lines for L_{d} = 12.71 and L_{d} = 21.00 in the experimental data. At low Q values (Ca ≤ 2.5 × 10^{−5}), there was little difference in splitting ratio for L_{d} ≈ 12 compared with L_{d} ≈ 20 (see Fig. 7). The theoretical results show that increasing L_{d} decreases L_{b}/L_{t} for all Ca. The theory predicts the existence of Ca_{crit} for all L_{d} and indicates that Ca_{crit} increases as L_{d} increases.
DISCUSSION
The results demonstrate that a long bubble passing through an individual bifurcation splits more evenly as flow, indicated by Ca, increases. For slowenough flows, e.g., Ca < Ca_{crit}, the bubble does not split, and the entire bubble enters the upper branch. Higher θ and higher Bo* lead to less even splitting. The bubble splitting behavior is complicated and depends on the effects of surface tension, inertia, viscosity, gravity, and pressure. The theoretical model represents a simplified version of the experimental situation but provides good qualitative agreement with the experimental results and explains the mechanisms involved in the observed behavior. Equation 13 indicates that the velocities of the bubble noses in the upper and lower daughter tubes sum to a constant, so if the upper one speeds up, the lower one slows down. Higher daughter tube velocities and larger values of l result in higher viscous losses, as indicated by Eqs. 3 and 4. Equations 12 and 13 indicate that dl_{1}/dt and dl_{2}/dt, e.g., −V_{1} and −V_{2}, respectively, each have terms that depend on flow and a term that depends on gravity. When θ = 0°, the gravitational term is zero and V_{1} = V_{2}, indicating even bubble splitting. When θ > 0°, the gravitational term is nonzero and V_{1} > V_{2}, indicating that more of the bubble enters the upper daughter tube. However, the effect of Q increases with Q, and for very high Q values V_{1} ≈ V_{2}, indicating nearly even splitting. For slowenough flow or highenough θ, the bubble cannot overcome gravity's opposition to its motion into the lower daughter tube, and the splitting ratio is zero. Higher values of L_{d} correspond to longer final times in integrating Eqs. 12 and 13, resulting in larger changes in l_{1} and smaller splitting ratios, as in the experiments.
This quasisteady analysis captures the competition between flow and gravity, but neglects the specific details of the interface shape and the threedimensional features of the flow. Consequently, it does not exhibit the stickslip behavior observed in the arteriole model. The theoretical model does, however, capture the bubble reversal phenomenon that was observed in the experiments. As indicated by Eq. 2, gravity opposes the motion of the bubble into the lower branch. For high Q values, the gravitational effects are negligible compared with the viscous pressure drop along the tube. For lower Q values, the effects of gravity are significant, and, as the menisci in the upper and lower daughter branches propagate, the elevation difference between them grows, causing the upper meniscus to accelerate and the lower one to decelerate. One can see this from Eqs. 3, 4, 8, and 9, by considering the simplified system in which the velocity dependence of P_{1} and P_{2} in Eq. 7 is neglected and l_{1} and l_{2} are negligible compared with L. In that case, P_{1} = P_{2}, and an increase in the elevation difference between points 1 and 2, which occurs as the bubble moves into the bifurcation, is accompanied by an increase in the difference between V_{1} and V_{2}. Mass conservation, Eq. 10, requires V_{1} and V_{2} to sum to a constant, so that, as the upper meniscus accelerates, the lower one decelerates. Once V_{2} becomes negative, the lower meniscus moves backward out of the lower daughter tube and toward the parent tube. If there is sufficient time before the bubble splits, this can result in L̂_{b}(t̂) reaching zero and the bubble moving entirely into the upper daughter branch. Without those simplifications, the behavior is more complicated, but bubble reversal occurs as shown in Fig. 5. The theoretical model results agree well with the artery model results. Qualitatively, the theory predicts the arteriole model experimental behavior, but quantitative agreement is not as good. Due to the smaller Bo* in the arteriole model, the effects of surface tension are more important in the arteriole model than in the artery model. A more complicated, computational model that accurately predicts the bubble interface shape and motion is needed to quantitatively capture the behavior in the arteriole model.
These findings explain the air emboli behavior observed by Chang et al. (12) and Souders et al. (41). In vessels in which Ca is less than Ca_{crit} and Bo* is high, such as the pulmonary artery and the first several subsequent vessel generations, bubbles will not divide and will be transported to the upper branches. This results in the transport of air emboli to the superior regions of the lungs, where they can lodge, as demonstrated by Chang et al. (12). However, when emboli reach smaller vessels, we expect them to split because Ca will be greater than Ca_{crit} for the particular value of Bo* in small vessels, which scales as vessel diameter squared and is small. The splitting behavior in these arterioles corresponds to the more even splitting observed in our arteriole model. In even smaller vessels, corresponding to yet smaller Bo*, the emboli will split more evenly. This suggests that, in small pulmonary vessels, the importance of flow in determining the emboli distribution will be greater than in large vessels. This explains, at the size scale of a single bifurcation, why Souders et al. (41) observed that bubbles introduced into the lungs of anesthetized dogs were influenced by blood flow as well as gravity. Although they did not consider bubble transport at the bifurcation level, they concluded that flow was a factor in the distribution of emboli within the lung, based on their measurements of regional blood flow distribution. The findings of Souders et al. and Chang et al. (12) are consistent in that they both demonstrated the propensity of emboli to enter the superior regions of the lungs. However, Souders et al. (41) observed that emboli transport in small vessels was not entirely determined by gravity. The present study demonstrates that gravity may dominate in large vessels, such as those resolved by the xenon tracer method used by Chang et al. (12), but that flow can dominate in smaller vessels, which could be resolved by the more accurate fluorescent microsphere method used by Souders et al. (41).
In a related study, Chang et al. (13) investigated the transport of small (L_{d} << 1) bubbles in a bifurcation with a 90° θ. Their experiments suggested that a small bubble will enter the upper or lower branch, depending on its vertical position in the parent tube before reaching the bifurcation, and they noted that the vast majority of the bubbles went to the upper branch because they were near the top of the parent tube before entering the bifurcation. Based on this, Chang et al. (13) concluded that Q did not significantly influence the behavior of the air emboli. Manga (37) computationally investigated the motion of small (L_{d} << 1) liquid droplets transported by a suspending liquid through a bifurcation without gravity as a model of geological flows in porous media. The value of Q was specified at the inlet and both outlets of the bifurcation. The small droplets did not split and generally moved into the daughter branch with the higher Q (37). In vivo studies have indicated that air emboli in microvessels are typically not small compared with vessel diameter, rather they are long, sausageshaped bubbles (4). In the present study of long bubbles, splitting favored the upper branch at low Q values. Here, it was demonstrated that increasing Q resulted in more even splitting and that Q significantly affects the splitting ratio.
Due to the small size of capillaries, Bo* is much smaller in the microcirculation than in arteries, and, consequently, the value of Ca_{crit} is much smaller there. This suggests that, if PFC bubbles for gas embolotherapy are formed in the microcirculation, they will split more evenly at each bifurcation than they would in arteries, as indicated by comparing the artery and arteriole model results, leading to a relatively homogenous distribution of PFC emboli. The choice of bubble size relative to the vessel diameter will likely affect homogeneity, as indicated by the bubble length, which is approximately six times the parent tube diameter, leading to more homogenous splitting than the longer bubbles in this study. It is not clear at this time if the potential for vessel rupture during vaporization will require formation of PFC bubbles in vessels larger than capillaries (44). If bubbles are formed in larger vessels, gravity will likely play an important role in determining their transport until they reach smaller vessels.
Although this idealized experimental model captures the effects of flow, gravity, surface tension, and viscosity, it has some limitations. For example, actual vessels are flexible and are lined by an endothelial layer, which interacts with air bubbles to change adhesion strength (42, 43). Blood is nonNewtonian and contains cellular components and surfaceactive proteins and lipids. Shear thinning of blood likely modifies the velocity and stress fields, both of which influence bubble transport, compared with the constant viscosity (Newtonian) fluid that we used to model blood. Surfactants have been shown to modify bubble sticking (19), to attenuate gas embolisminduced thrombin production (21), and to stabilize microbubbles (15, 16, 45). Computational studies of semiinfinite bubbles in rigid tubes have demonstrated that surfactant properties can influence the interfacial pressure drop through modification of the surface tension and through Marangoni stresses that affect the viscous stresses along the interface (26). Surfactant modifies the thickness of the liquid layer between the bubble and the tube wall (26). Surfactants have also been shown to have similar effects on finitelength bubbles in inclined tubes (17) and to modify the dynamics of liquid droplets (22–24, 39). Surfaceactive proteins and lipids in the blood likely affect bubble splitting by similar mechanisms.
Flow in arteries is pulsatile, but only steady flow was considered here. Our experimental setup did not allow Q values higher than 80 ml/min, at which the bubble moved too fast for accurate imaging via the camera. Another limitation was the inability to effectively vary the bubble size in the artery model, due to bubble break up. Therefore, we only had one bubble size for the artery model. Although it is unlikely that larger bubbles or larger arteries would be targeted for gas embolotherapy, they may be of interest with regard to air embolism. Future studies may incorporate some of the physiological features that have been neglected in this investigation and should consider the effects of bifurcation asymmetry, either from geometry or a pressure difference between the daughter branches. Tumor microcirculation is complicated and irregular, and, consequently, more in vitro and in vivo experiments are needed to understand bubble transport there. This study considered constant Q values, and, therefore, did not consider bubble sticking, which is not expected to occur in the vessels modeled here (40). Previous computational (31) and experimental (27) studies have suggested that the pressure jump across a bubble interface can be affected by inertia but do not provide an analytical, theoretically based relationship for the inertial contribution to ΔP. Inertial effects have been neglected in the present model due to the dominance of the viscous pressure drop in the unidirectional Poiseuille flow ahead of the bubble, and previous computational studies (28) suggest that the effects of inertia on the film thickness of a semiinfinite bubble are minimal for Ca < 0.05. However, these effects could be included in future models to allow investigation of a wider parameter range with improved accuracy. Timedependent computational models of bubble transport in bifurcations are planned to investigate the dynamics, such as the bubble snapoff and velocity fields near the bubble nose, that were not captured by the quasisteady model presented here.
In conclusion, it was shown that the splitting ratio increases with capillary number (flow) and decreases with θ. A Ca_{crit} exists, below which the bubble does not divide and entirely enters the upper branch. The value of the Ca_{crit} increases with the effective Bo. This interplay of gravity and flow at the bifurcation level is likely responsible for the previously observed distribution of gas emboli in the lungs, in which emboli typically lodge in the superior regions of the lungs but have a flowdependent distribution (12, 41) within smaller vessels. For the long bubbles considered in our arteriole model, those with L_{d} ≈ 6 split more evenly than longer bubbles (L_{d} ≈ 12 and 20). The findings of this study may potentially be incorporated into an embolotherapy strategy to achieve a relatively homogenous occlusion of blood flow to tumors.
GRANTS
This work is funded by National Science Foundation Grant BES0301278 and National Institutes of Health Grants EB003541 and EB00281.
Acknowledgments
We acknowledge the assistance of John Myers in the construction of the bifurcations and helpful discussions with Dr. Oliver Kripfgans and Catherine Orifici.
Footnotes
The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “advertisement” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
 Copyright © 2005 the American Physiological Society